﻿namespace EnergyGeneratorSim
{
    using System;
    using System.IO;
    using EnergyGeneratorSim.Parameters;

    public class RaumenergiemotorDFEM
    {
        const int Bn = 7;          // Schritt-Anzahl der Magnetfeld-Speicherung nach *2 von S.2
        const int SpNmax = 200;    // Maximal mögliche Anzahl der Stützpunkte der Spulen (Inupt und Turbo)
        const int FlNmax = 2000;   // Maximal mögliche Anzahl der Flächenelemente der Spulen (Inupt und Turbo)
        const int MESEanz = 200;   // Tatsächliche Anzahl der Magnetfeld-Emulations-Spulen-Elemente, gerade Anzahl wählen!
        const int AnzPmax = 35000; // Dimensionierung der Arrays für den Plot (für den Datenexport nach Excel)

        //private double Bsw;  // Schritt-Weite der Magnetfeld-Speicherung nach *1 von S.2
        private double Spsw;  // Schritt-Weite der Spulen-Aufgliederung nach *2 von S.1
        private int SpN; // Anzahl der Stützpunkte der Spulen
        private int FlN; // Anzahl der Flächenelemente der Spulen
        private int xo, yo, zo; // Geometrieparameter nach Zeichnung *2 von S.1
        private double PsiSFE;  // magnetischer Fluß durch ein Spulen-Flächenelement
        private double PsiGES;  // magnetischer Fluß durch die gesamte Spule
        private double B1T, B2T, B3T, B4T, B5T;  // Fourier-Koeffizienten, Turbo-Spule
        private double B1I, B2I, B3I, B4I, B5I;  // Fourier-Koeffizienten, Input-Spule
        private double B1dreh, phase;  // Koeffizienten zur Drehmoments-Schnellberechnung
        //private double MEyo, MEro, MEI;// Abmessungen und Strom des Magnetfeld-Emulationsspulenpaares nach *1 von S.5
        private double[, ,] Bx = new double[2 * Bn + 1, 2 * Bn + 1, 2 * Bn + 1];// Kartes. Komp. der magn. Induktion, Dauermagnet nach *3 von S.2
        private double[, ,] By = new double[2 * Bn + 1, 2 * Bn + 1, 2 * Bn + 1];
        private double[, ,] Bz = new double[2 * Bn + 1, 2 * Bn + 1, 2 * Bn + 1];
        private double[, ,] OrtBx = new double[2 * Bn + 1, 2 * Bn + 1, 2 * Bn + 1];// Kartes. Komp. der Orte, an denen das Feld Bx, By, Bz ist.
        private double[, ,] OrtBy = new double[2 * Bn + 1, 2 * Bn + 1, 2 * Bn + 1];
        private double[, ,] OrtBz = new double[2 * Bn + 1, 2 * Bn + 1, 2 * Bn + 1];
        private double[, ,] BxDR = new double[2 * Bn + 1, 2 * Bn + 1, 2 * Bn + 1]; // gedrehtes Magnetfeld
        private double[, ,] OrtBxDR = new double[2 * Bn + 1, 2 * Bn + 1, 2 * Bn + 1];
        private double[, ,] OrtByDR = new double[2 * Bn + 1, 2 * Bn + 1, 2 * Bn + 1];
        private double[, ,] OrtBzDR = new double[2 * Bn + 1, 2 * Bn + 1, 2 * Bn + 1]; // gedrehte Ortsvektoren
        private double[] MESEx = new double[MESEanz];// Orte der Magnetfeld-Emulations-Spulen-Elemente
        private double[] MESEy = new double[MESEanz];
        private double[] MESEz = new double[MESEanz];
        private double[] MESEdx = new double[MESEanz];// Laufrichtungen der Magnetfeld-Emulations-Spulen-Elemente
        private double[] MESEdy = new double[MESEanz];
        private double[] MESEdz = new double[MESEanz];
        private double[] SpIx = new double[SpNmax - 1];// Stützpunkte der Polygonzüge der Input-Spule, kartesische Koordinaten
        private double[] SpIy = new double[SpNmax - 1];
        private double[] SpIz = new double[SpNmax - 1];
        private double[] SpTx = new double[SpNmax - 1];// Stützpunkte der Polygonzüge der Turbo-Spule, kartesische Koordinaten
        private double[] SpTy = new double[SpNmax - 1];
        private double[] SpTz = new double[SpNmax - 1];
        private double[] SIx = new double[SpNmax - 1];// Ort als Mittelpunkt der Leiterschleifen-Elemente
        private double[] SIy = new double[SpNmax - 1];
        private double[] SIz = new double[SpNmax - 1];
        private double[] STx = new double[SpNmax - 1];// Ort als Mittelpunkt der Leiterschleifen-Elemente
        private double[] STy = new double[SpNmax - 1];
        private double[] STz = new double[SpNmax - 1];
        private double[] dSIx = new double[SpNmax - 1];// Richtungsvektoren der Leiterschleifen-Elemente der Input-Spule
        private double[] dSIy = new double[SpNmax - 1];
        private double[] dSIz = new double[SpNmax - 1];
        private double[] dSTx = new double[SpNmax - 1];// Richtungsvektoren der Leiterschleifen-Elemente der Turbo-Spule
        private double[] dSTy = new double[SpNmax - 1];
        private double[] dSTz = new double[SpNmax - 1];
        private double[] FlIx = new double[SpNmax - 1];// Flächenelemente der Input-Spule, kartesische Koordinaten
        private double[] FlIy = new double[SpNmax - 1];
        private double[] FlIz = new double[SpNmax - 1];
        private double[] FlTx = new double[SpNmax - 1];// Flächenelemente der Turbo-Spule, kartesische Koordinaten
        private double[] FlTy = new double[SpNmax - 1];
        private double[] FlTz = new double[SpNmax - 1];
        // Zum Lösen der Bewegungs-Differentialgleichung:
        private double phio, phipo, phippo, phim, phipm; // Winkel und dessen Ableitungen Index "io" und "io-1"
        private double qoT, qpoT, qppoT, qmT, qpmT; // Ladung und deren Ableitungen Index "io" und "io-1" in Turbo-Spule
        private double qoI, qpoI, qppoI, qmI, qpmI, qppmI; // Ladung und deren Ableitungen Index "io" und "io-1" in Input-Spule
        private double UinduzT, UinduzI;   // Induzierte Spannung zum "Jetzt-Zeitpunkt" in Input-Spule und Turbo-Spule
        private ulong znr;  // Zählnummer für's Plotten der Daten ins Excel.
        private ulong LPP; // der Wert wird für Datenausgabe benutzt.

        private SimulationDataItem[] simPrep = new SimulationDataItem[360];

        private double BTx, BTy, BTz;  // Magnetfeld der Turbo-Spule an einem beliebigen Aufpunkt
        private double BIx, BIy, BIz;  // Magnetfeld der Input-Spule an einem beliebigen Aufpunkt
        private double merk;  // Für Test-Zwecke
        private double DLI, DLT;  // Durchmesser und Länge des Spulendrahtes zur Angabe der Drahtstärke
        private double RI, RT;  // Ohm'scher Widerstand der Spulen-Drähte
        private double LT, LI;  // Induktivität der Turbo-Spule und der Input-Spule
        private double BreiteI, HoeheI, BreiteT, HoeheT; // Breite und Höhe der beiden Spulenlörper
        private double omT, TT;  // Kreis-Eigenfrequenz und Schwingungsdauer des Turbo-Spulen-Schwingkreises aus LT & CT.
        private double omAn;  // Umdrehungen pro Minute Startdrehzahl, und Winkelgeschwindigkeit rad/sec. Startdrehzahl
        private double UmSec;  // Umdrehungen pro Sekunde Startdrehzahl
        private double J;  // Trägheitsmoment des Magneten bei Rotation
        private double Mmag;  // Masse des Magneten
        private double Tjetzt;  // Aktueller Zeitpunkt (beim Lösen der Dgl.)
        private double QTmax, QImax, QpTmax, QpImax, QppTmax, QppImax, phipomax; // Maximalwerte finden für Strom-, Drehzahl- und Spannungsangabe
        private double Wentnommen; // entnommene Gesamtenergie
        private double AnfEnergie, EndEnergie; // Energie-Vergleich im System
        private bool steigtM, steigtO;  // steigt die Flanke eines Referenz-Signals (als und neu) ?
        private ulong Iumk;  // Lage des Umkehrpunktes als Referenz für das Input-Spannungs-Signal
        private double fkI, fkT; // Korrekturfaktoren zur Induktivität lt. Stöcker S. 452
        private double Pzuf, Ezuf; // Zugeführte Leistung und Energie über die Input-Spannung
        private double cr; // Reibungs-Koeffizient für geschwindigkeits-proportionale Reibung
        private double Preib; // Über mechanische Reibung entnommene Leistung
        private double Ereib; // Entnommene mechanische Energie, z.B. über Reibung

        /// <summary>
        /// Initializes a new instance of the <see cref="RaumenergiemotorDFEM"/> class.
        /// </summary>
        public RaumenergiemotorDFEM()
        {
            for (int i = 0; i < simPrep.Length; i++)
            {
                simPrep[i] = new SimulationDataItem();
            }

            this.ProgressReportStepCount = 1;
        }

        /// <summary>
        /// Gets or sets a value indicating whether this <see cref="RaumenergiemotorDFEM"/> is cancel.
        /// Set this to true to cancel the current calculation.
        /// </summary>
        /// <value><c>true</c> if cancel; otherwise, <c>false</c>.</value>
        public bool Cancel { get; set; }

        private int MapArray(int i)
        {
            return MapArray(i, -Bn, Bn);
        }

        /// <summary>
        /// Maps the array index from the range lowerLimit..upperLimit to the 
        /// range 0 .. (lowerLimit+upperLimit)+1.
        /// </summary>
        /// <param name="i">The loop index.</param>
        /// <param name="lowerLimit">The lower limit.</param>
        /// <param name="upperLimit">The upper limit.</param>
        /// <returns></returns>
        private static int MapArray(int i, int lowerLimit, int upperLimit)
        {
            if (i < lowerLimit || i > upperLimit)
            {
                throw new ArgumentOutOfRangeException("i", i, "Array index out of range!");
            }

            return i - lowerLimit;
        }

/*
        private void Dokumentation_des_Ergebnisses(SimulationParameters param)
        {
            using (var stream = File.OpenWrite("Auswertung"))
            using (var fout = new StreamWriter(stream))
            {
                fout.WriteLine("DFEM-Simulation eines EMDR-Motors.");
                fout.WriteLine(" ");
                fout.WriteLine("Parameter zum Lösen der Dgl. und zur Darstellung der Ergebnisse:");
                fout.WriteLine("AnzP = " + param.AnzP + " Zum Lösen der Dgl.: Anzahl der tatsächlich berechneten Zeit-Schritte");
                fout.WriteLine("dt   = " + param.dt + " // Sekunden Dauer der Zeitschritte zur iterativen Lsg. der Dgl.");
                fout.WriteLine("Abstd= " + param.Abstd + "  {Nur Vorbereitung, nicht zum Lösen der Dgl.: Jeden wievielten Punkt soll ich plotten ins Excel?");
                fout.WriteLine("PlotEnde   = " + param.PlotEnde + " {Zum Lösen der Dgl.: Letzter-Plot-Punkt: Ende des Daten-Exports nach Excel");
                fout.WriteLine("PlotStep   = " + param.PlotStep + " {Zum Lösen der Dgl.: Schrittweite des Daten-Exports nach Excel");
                fout.WriteLine(" ");
                fout.WriteLine("PlotAnfang = " + param.PlotAnfang + " {Zum Lösen der Dgl.: Erster-Plot-Punkt: Anfang des Daten-Exports nach Excel");
                fout.WriteLine("Es folgt die Eingabe der beiden Spulen, vgl. Zeichnung *2 von S.1 :");
                fout.WriteLine("Die Spulen werden nach Vorgabe der Geometrieparameter automatisch vernetzt.");
                fout.WriteLine("Spsw   = " + Spsw + " Angabe in Metern: Die Spulen-Aufgliederung ist in Spsw-Schritten");
                fout.WriteLine("xo = " + xo + ",  Angaben in Vielfachen von Spsw, Geometrieparameter nach Zeichnung*2 von S.1");
                fout.WriteLine("yo = " + yo + ",  Angaben in Vielfachen von Spsw, Geometrieparameter nach Zeichnung*2 von S.1");
                fout.WriteLine("zo = " + zo + ",  Angaben in Vielfachen von Spsw, Geometrieparameter nach Zeichnung*2 von S.1");
                fout.WriteLine("Ninput = " + param.Inductor.Ninput + "  Zahl der Wicklungen der Input-Spule");
                fout.WriteLine("Nturbo = " + param.Inductor.Nturbo + "  Zahl der Wicklungen der Turbo-Spule");
                fout.WriteLine("nebeninput = " + param.Inductor.Nebeninput + "  Windungen nebeneinander in der Input-Spule");
                fout.WriteLine("ueberinput = " + param.Inductor.Ueberinput + "  Windungen übereinander in der Input-Spule");
                fout.WriteLine("nebenturbo = " + param.Inductor.Nebenturbo + "  Windungen nebeneinander in der Turbo-Spule");
                fout.WriteLine("ueberturbo = " + param.Inductor.Ueberturbo + "  Windungen übereinander in der Turbo-Spule");
                fout.WriteLine(" ");
                fout.WriteLine("Bsw = " + param.PermanentMagnet.Bsw + "  Magnetfeld-Speicherung nach *1 von S.2 in Zentimeter-Schritten");
                fout.WriteLine("Ich emuliere hier das Magnetfeld eines 1T-Magneten durch ein Spulenpaar nach *1 von S.5:");
                fout.WriteLine("MEyo = " + param.PermanentMagnet.MEyo + "  y-Koordinaten der Magnetfeld-Emulationsspulen nach *1 von S.5");
                fout.WriteLine("MEro = " + param.PermanentMagnet.MEro + "  Radius der Magnetfeld-Emulationsspulen nach *1 von S.5");
                fout.WriteLine("MEI  = " + param.PermanentMagnet.MEI + "  Strom des Magnetfeld-Emulationsspulenpaares nach *1 von S.5, Angabe in Ampere");
                fout.WriteLine(" ");
                fout.WriteLine("Allgemeine technische Größen:");
                fout.WriteLine("DD = " + param.Inductor.DD + "  // Meter {Durchmesser des Spulendrahtes zur Angabe der Drahtstärke");
                fout.WriteLine("rho = " + param.General.rho + "  // Ohm*m    {Spez. elektr. Widerstand von Kupfer, je nach Temperatur, Kohlrausch,T193");
                fout.WriteLine("rhoMag = " + param.General.rhoMag + "  // kg/m^3 {Dichte des Magnet-Materials, Eisen, Kohlrausch Bd.3");
                fout.WriteLine("CT = " + param.General.CT + "  // Farad    {Kapazität des Kondensators, der mit in der Turbo-Spule (!) in Reihe geschaltet");
                fout.WriteLine("CI = " + param.General.CI + "  // Farad    {Kapazität des Kondensators, der mit in der Input-Spule (!) in Reihe geschaltet");
                fout.WriteLine(" ");
                fout.WriteLine("Sonstige (zur Eingabe):");
                fout.WriteLine("Rlast = " + param.General.Rlast + "  // Ohm   Ohm´scher Lastwiderstand im LC-Turbo-Schwingkreis");
                fout.WriteLine("UmAn = " + param.General.UmAn + "  // U/min Anfangsbedingung mechanisch - Rotierender Magnet: Startdrehzahl");
                fout.WriteLine("Uc = " + param.General.Uc + "  // Volt Anfangsbedingung elektrisch - Kondensatorspannung am TURBO-Kondensator");
                fout.WriteLine("Il = " + param.General.Il + "  // Ampere Anfangsbedingung elektrisch - Spulenstrom im TURBO-Schwingkreis");
                fout.WriteLine(" ");
                fout.WriteLine("Mechanische Leistungs-Entnahme, proportional zur Geschwindikeit, aber mit Nachregelung zur konst. Drehzahl:");
                fout.WriteLine("Koeffizient einer geschw-prop. mechan. Leistungs-Entnahme: " + param.General.crAnfang + " Nm/(rad/s)");
                fout.WriteLine("Ziel-Drehzahl,für mechanische Reibungs-Nachregelung: " + param.General.phipZiel + " U/min");
                fout.WriteLine(" ");
                fout.WriteLine("Abgeleitete Parameter. Die Größen werden aus den obigen Parametern berechnet, es ist keine Eingabe möglich:");
                fout.WriteLine("DLI:=4*(yo+zo)*Spsw*Ninput = " + DLI + " // Meter Länge des Spulendrahtes, Input-Spule");
                fout.WriteLine("DLT:=4*(yo+zo)*Spsw*Nturbo = " + DLT + " // Meter Länge des Spulendrahtes, Turbo-Spule");
                fout.WriteLine("RI:=rho*(DLI)/(Math.PI/4*DD*DD) = " + RI + " // Ohm Ohm`scher Widerstand des Spulendrahtes, Input-Spule");
                fout.WriteLine("RT:=rho*(DLT)/(Math.PI/4*DD*DD) = " + RT + " // Ohm Ohm`scher Widerstand des Spulendrahtes, Turbo-Spule");
                fout.WriteLine("BreiteI:=nebeninput*DD = " + BreiteI + " Breite und Höhe des Input-Spulenlörpers");
                fout.WriteLine("HöheI:=ueberinput*DD  = " + HoeheI + " Breite und Höhe des Input-Spulenlörpers");
                fout.WriteLine("BreiteT:=nebenturbo*DD = " + BreiteT + " Breite und Höhe des Turbo-Spulenlörpers");
                fout.WriteLine("HöheT:=ueberturbo*DD  = " + HoeheT + " Breite und Höhe des Turbo-Spulenlörpers");
                fout.WriteLine("fkI:=Math.Sqrt(HöheI*HöheI+4/Math.PI*2*yo*2*zo)/HöheI = " + fkI + " Korrekturfaktor zur Induktivität der kurzen Input-Spule");
                fout.WriteLine("fkT:=Math.Sqrt(HöheT*HöheT+4/Math.PI*2*yo*2*zo)/HöheT = " + fkT + " Korrekturfaktor zur Induktivität der kurzen Turbo-Spule");
                fout.WriteLine("LI:=muo*(2*yo+BreiteI)*(2*zo+BreiteI)*Ninput*Ninput/(HöheI*fkI) = " + LI + " Induktivität Input-Spule");
                fout.WriteLine("LT:=muo*(2*yo+BreiteT)*(2*zo+Breitet)*Nturbo*Nturbo/(HöheT*fkT) = " + LT + " Induktivität Turbo-Spule");
                fout.WriteLine("omT:=1/Math.Sqrt(LT*CT) = " + omT + " Kreis-Eigenfrequenz des Turbo-Spulen-Schwingkreises aus LT & CT");
                fout.WriteLine("TT:=2*Math.PI/omT = ", TT, " Schwingungsdauer des Turbo-Spulen-Schwingkreises aus LT & CT.");
                fout.WriteLine("Mmag:=rhoMag*(Math.PI*MEro*MEro)*(2*MEyo) = " + Mmag + " kg Masse des Magneten // Rotation um Querachse !!");
                fout.WriteLine("J:=Mmag/4*(MEro*MEro+4*MEyo*MEyo/3) = " + J + " Trägheitsmoment des Magneten bei Rotation um Querachse");
                fout.WriteLine(" ");
                fout.WriteLine("Anzeige einiger auszurechnender Parameter:");
                fout.WriteLine("Magnet: Start-Winkelgeschw.: omAn = " + omAn + " rad/sec");
                fout.WriteLine("Magnet: Startdrehzahl, Umdr./sec.: UmSec = " + UmSec + " Hz");
                fout.WriteLine("Masse des Magnet = " + Mmag + " kg");
                fout.WriteLine("Tr'gheitsmoment Magnet bei QUER-Rotation" + J + " kg*m^2");
                fout.WriteLine("Gesamtdauer der Betrachtung: " + param.AnzP * param.dt + " sec.");
                fout.WriteLine("Excel-Export: " + param.PlotAnfang * param.dt + "..." + param.PlotEnde * param.dt + " sec., Step " + param.PlotStep * param.dt + " sec.");
                fout.WriteLine("Das sind " + (param.PlotEnde - param.PlotAnfang) / param.PlotStep + " Datenpunkte (also Zeilen).");
                fout.WriteLine(" ");
                fout.WriteLine(" *********************************************************************************************************");
                fout.WriteLine(" ");
                fout.WriteLine("Einige Ergebnisse der Berechnung:");
                fout.WriteLine("Anfangs-Energie im System:    " + AnfEnergie + " Joule");
                fout.WriteLine("End-Energie im System:        " + EndEnergie + " Joule");
                fout.WriteLine("Leistungs-Änderung im System:" + (EndEnergie - AnfEnergie) / (param.AnzP * param.dt) + " Watt");
                fout.WriteLine("Am Lastwiderstand entnommene Gesamtenergie = " + Wentnommen + " Joule");
                fout.WriteLine("entsprechend einer mittleren entnommenen Leistg:" + Wentnommen / (param.AnzP * param.dt) + " Watt");
                fout.WriteLine("Über Input-Spannung zugeführte Gesamt-Energie: " + Ezuf + " Joule");
                fout.WriteLine("entsprechend einer mittleren zugeführten Leistg:" + Ezuf / (param.AnzP * param.dt) + " Watt");
                fout.WriteLine("Gesamte mechanisch entnommene Energie = " + Ereib + " Joule");
                fout.WriteLine("entsprechend einer mittleren Leistung = " + Ereib / (param.AnzP * param.dt) + " Watt");
                fout.WriteLine("bei einer Betrachtungs-Dauer von " + (param.AnzP * param.dt) + " sec.");
            }
        }
*/

        private void ExcelAusgabe(SimulationParameters param, SimulationDataItem[] simData, string Name, int Spalten)
        {
            using (var stream = File.OpenWrite(Name))
            using (var fout = new StreamWriter(stream))
            {
                for (ulong lv = 0; lv < param.AnzP; lv++)
                {
                    if (lv % param.Abstd != 0)
                        continue;

                    for (int j = 0; j < Spalten; j++)
                    {
                        // Kolumnen drucken, zuerst 3*Ladung, dann 3*Winkel, dann 8 freie Felder
                        double zahl = double.NaN;
                        if (j == 0) zahl = simData[lv].Q;
                        if (j == 1) zahl = simData[lv].Qp;
                        if (j == 2) zahl = simData[lv].Qpp;
                        if (j == 3) zahl = simData[lv].phi;
                        if (j == 4) zahl = simData[lv].phip;
                        if (j == 5) zahl = simData[lv].phipp;
                        if (j == 6) zahl = simData[lv].KG;
                        if (j == 7) zahl = simData[lv].KH;
                        if (j == 8) zahl = simData[lv].KI;
                        if (j == 9) zahl = simData[lv].KJ;
                        if (j == 10) zahl = simData[lv].MagneticFluxTurbo;
                        if (j == 11) zahl = simData[lv].MageneticFluxInput;
                        if (j == 12) zahl = simData[lv].VoltageInducedTurbo;
                        if (j == 13) zahl = simData[lv].VoltageInducedInput;

                        fout.Write(zahl.ToString().Replace(".", ","));
                        fout.Write("\t");  // Daten-Trennung, Tabulator
                    }

                    fout.WriteLine();    // Zeilen-Trennung
                }
            }
        }

        private void ExcelLangAusgabe(SimulationParameters param, SimulationDataItem[] simData, string Name, int Spalten)
        {
            if (Spalten > 25)
            {

                PrintLine("FEHLER: Zu viele Spalten. Soviele Daten-Arrays sind nicht vorhanden.");
                PrintLine(" => PROGRAMM WURDE ANGEHALTEN : STOP !");
                Environment.Exit(1);
            }

            using (var stream = File.OpenWrite(Name))
            using (var fout = new StreamWriter(stream))
            {
                // Überschriften andrucken
                fout.WriteLine("Zeit\tQ\tQp\tQpp\tQI\tQpI\tQppI\tphi\tphip\tphipp\tFlussT\tFlussI\tUindT\tUindI\tEinputSp\tEturboSp\tEinputKond\tEturboKond\tEges\tElast\tUinput\tEzu\tcr\tEmech\tleer");

                for (ulong lv = 0; lv < LPP; lv++)   // von "plotanf" bis "plotend"
                {
                    if ((lv % param.Abstd) != 0)
                        continue;

                    for (int j = 0; j < Spalten; j++)
                    {   // Kolumnen drucken, zuerst 3*Ladung, dann 3*Winkel, dann 8 freie Felder
                        double zahl = double.NaN;

                        if (j == 0) zahl = simData[lv].Zeit;  // Markieren der Zeit-Skala
                        if (j == 1) zahl = simData[lv].Q;     // Turbo-Spule
                        if (j == 2) zahl = simData[lv].Qp;    // Turbo-Spule
                        if (j == 3) zahl = simData[lv].Qpp;   // Turbo-Spule
                        if (j == 4) zahl = simData[lv].QI;    // Input-Spule
                        if (j == 5) zahl = simData[lv].QpI;   // Input-Spule
                        if (j == 6) zahl = simData[lv].QppI;  // Input-Spule
                        if (j == 7) zahl = simData[lv].phi;   // Magnet
                        if (j == 8) zahl = simData[lv].phip;  // Magnet
                        if (j == 9) zahl = simData[lv].phipp; // Magnet
                        if (j == 10) zahl = simData[lv].MagneticFluxTurbo;    // Auxiliary
                        if (j == 11) zahl = simData[lv].MageneticFluxInput;    // Auxiliary
                        if (j == 12) zahl = simData[lv].VoltageInducedTurbo;    // Auxiliary
                        if (j == 13) zahl = simData[lv].VoltageInducedInput;    // Auxiliary
                        if (j == 14) zahl = simData[lv].EnergyInputCoil;    // Auxiliary
                        if (j == 15) zahl = simData[lv].EnergyTurboCoil;    // Auxiliary
                        if (j == 16) zahl = simData[lv].EnergyInputCapacitor;    // Auxiliary
                        if (j == 17) zahl = simData[lv].EnergyTurboCapacitor;    // Auxiliary
                        if (j == 18) zahl = simData[lv].EnergyMagneticRotation;    // Auxiliary
                        if (j == 19) zahl = simData[lv].EnergySystem;    // Auxiliary
                        if (j == 20) zahl = simData[lv].PowerElectric;    // Auxiliary
                        if (j == 21) zahl = simData[lv].VoltageInputCoil;    // Auxiliary
                        if (j == 22) zahl = simData[lv].PowerAddedByInput;    // Auxiliary
                        if (j == 23) zahl = simData[lv].FrictionCoefficient;    // Auxiliary
                        if (j == 24) zahl = simData[lv].MechanicalFrictionPowerOut;    // Auxiliary
                        if (j == 25) zahl = 0;    // Auxiliary

                        fout.Write(zahl.ToString().Replace(".", ","));
                        fout.Write("\t");  // Daten-Trennung, Tabulator
                    }

                    fout.WriteLine();    // Zeilen-Trennung
                }
            }
        }

        /*
        private void Magnetfeld_zuweisen_01(PermanentMagnetParameters permanentMagnet) // homogenes Magnetfeld
        {
            for (int i = -Bn; i < Bn; i++)  // in x-Richtung
            {
                for (int j = -Bn; j < Bn; j++)  // in y-Richtung
                {
                    for (int k = -Bn; k < Bn; k++)  // in z-Richtung
                    {
                        Bx[MapArray(i), MapArray(j), MapArray(k)] = 0.0; // Telsa
                        By[MapArray(i), MapArray(j), MapArray(k)] = 1.0; // Telsa
                        Bz[MapArray(i), MapArray(j), MapArray(k)] = 0.0; // Telsa
                        OrtBx[MapArray(i), MapArray(j), MapArray(k)] = i * permanentMagnet.Bsw;
                        OrtBy[MapArray(i), MapArray(j), MapArray(k)] = j * permanentMagnet.Bsw;
                        OrtBz[MapArray(i), MapArray(j), MapArray(k)] = k * permanentMagnet.Bsw;
                    }
                }
            }
        }
*/

/*
        private void Magnetfeld_zuweisen_02(PermanentMagnetParameters permanentMagnet) // willkürlicher Versuch eines inhomogenen Magnetfeldes
        {
            for (int i = -Bn; i < Bn; i++)  // in x-Richtung
            {
                for (int j = -Bn; j < Bn; j++)  // in y-Richtung
                {
                    for (int k = -Bn; k < Bn; k++)  // in z-Richtung
                    {
                        Bx[MapArray(i), MapArray(j), MapArray(k)] = -Sgn(i) / (i * i + j * j + k * k + 1);

                        if (i == 0)
                        {
                            Bx[MapArray(i), MapArray(j), MapArray(k)] = 0; // Telsa
                        }

                        By[MapArray(i), MapArray(j), MapArray(k)] = 10 / (i * i + j * j + k * k + 1);                           // Telsa
                        Bz[MapArray(i), MapArray(j), MapArray(k)] = -Sgn(k) / (i * i + j * j + k * k + 1);

                        if (k == 0)
                        {
                            Bz[MapArray(i), MapArray(j), MapArray(k)] = 0; // Telsa
                        }

                        OrtBx[MapArray(i), MapArray(j), MapArray(k)] = i * permanentMagnet.Bsw;
                        OrtBy[MapArray(i), MapArray(j), MapArray(k)] = j * permanentMagnet.Bsw;
                        OrtBz[MapArray(i), MapArray(j), MapArray(k)] = k * permanentMagnet.Bsw;
                    }
                }
            }
        }
*/

        private void Berechne_dH(PermanentMagnetParameters permanentMagnet, double om, double t, double sx, double sy, double sz, double qwill, double dphi, out double dHx, out double dHy, out double dHz)
        {
            double KRPx = -om * permanentMagnet.MEro * Math.Cos(om * t) * (permanentMagnet.MEyo - sy);
            double KRPy = +om * permanentMagnet.MEro * Math.Cos(om * t) * (permanentMagnet.MEro * Math.Cos(om * t) - sx) + om * permanentMagnet.MEro * Math.Sin(om * t) * (permanentMagnet.MEro * Math.Sin(om * t) - sz);
            double KRPz = -om * permanentMagnet.MEro * Math.Sin(om * t) * (permanentMagnet.MEyo - sy);
            double lmsbetrag = Physics.Sqr(permanentMagnet.MEro * Math.Cos(om * t) - sx) + Physics.Sqr(permanentMagnet.MEyo - sy) + Physics.Sqr(permanentMagnet.MEro * Math.Sin(om * t) - sz);
            lmsbetrag = Math.Sqrt(lmsbetrag);
            double lmsbetraghoch3 = lmsbetrag * lmsbetrag * lmsbetrag;

            if (lmsbetraghoch3 <= 1E-50)
            {
                dHx = 0;
                dHy = 0;
                dHz = 0;
            }
            else
            {
                dHx = qwill * KRPx / 4.0 / Math.PI / lmsbetraghoch3 * dphi / 2.0 / Math.PI;
                dHy = qwill * KRPy / 4.0 / Math.PI / lmsbetraghoch3 * dphi / 2.0 / Math.PI;
                dHz = qwill * KRPz / 4.0 / Math.PI / lmsbetraghoch3 * dphi / 2.0 / Math.PI;
            }
            //  PrintLine("Infinitesimales Feldelement: ",dHx:12:7,", ",dHy:12:7,", ",dHz:12:7," A/m"); 
        }

        private void BerechneHges(PermanentMagnetParameters permanentMagnet, out double Hx, out double Hy, out double Hz, double sx, double sy, double sz)
        {
            Hx = 0;
            Hy = 0;
            Hz = 0;  // Initialisierung des Gesamtfelds für die Addition der Feldelemente
            double qwill = 1;
            double om = 2.0 * Math.PI * permanentMagnet.MEI / qwill; // Ladung und der Kreisfrequenz in der Magnetfeld-Emulationsspule, *1 von S.7
            double dphi = 2.0 * Math.PI / 1000.0;    // Radianten bei der Aufteilung der Spule in 1000 Abschnitte

            for (int ilok = 0; ilok < 1000; ilok++)  // 1000 Zählschritte
            {
                double dHx, dHy, dHz;
                double t = ilok * dphi / om;  // Laufvariable (Zeit), zur Umrundung der Spule
                //    PrintLine("ilok = ",ilok:4," => ",om*t:12:6);  Wait; 
                Berechne_dH(permanentMagnet, om, t, sx, sy, sz, qwill, dphi, out dHx, out dHy, out dHz);   // Infinitesimales Feldelement nach Biot-Savart berechnen
                Hx = Hx + dHx;
                Hy = Hy + dHy;
                Hz = Hz + dHz;
            }

            //  PrintLine("Gesamtes Feld am Aufpunkt. : ",Hx:12:7,", ",Hy:12:7,", ",Hz:12:7," A/m"); 
            double Nenner;
            Nenner = Math.Sqrt(permanentMagnet.MEro * permanentMagnet.MEro + (permanentMagnet.MEyo - sy) * (permanentMagnet.MEyo - sy));
            Nenner = 2 * Nenner * Nenner * Nenner;
        }

        private void MagnetfeldZuweisen03(PermanentMagnetParameters permanentMagnet)
        {
            double sx, sy, sz; // Aufpunkt, an dem das Feld bestimmt werden soll
            double Hx, Hy, Hz; // Gesamtfeld am Ort des Aufpunkts
            double BXmax, BYmax, BZmax; // maximaler Feldwert auf der y-Achse

            PrintLine();
            PrintLine("Magnetfeld Emulations-Spulenpaar nach *1 von S.5");
            PrintLine("y-Koordinaten der Magnetfeld-Emulationsspulen nach *1 von S.5: " + permanentMagnet.MEyo + " m");
            PrintLine("Radius der Magnetfeld-Emulationsspulen nach *1 von S.5: " + permanentMagnet.MEro + " m");
            PrintLine("Strom der Magnetfeld-Emulationsspulen nach *1 von S.5: " + permanentMagnet.MEI + " Ampere");
            PrintLine("Anzahl der Schritte: " + Bn + " hoch 3 => " + 2 * Bn + 1 + " Bildschirm-Aktionspunkte je Spule.");

            //  Zuerst die obere Spule durchlaufen lassen: 
            for (int i2 = -Bn; i2 < Bn; i2++)  // in x-Richtung
            {
                for (int j2 = -Bn; j2 < Bn; j2++)  // in y-Richtung
                {
                    for (int k2 = -Bn; k2 < Bn; k2++)  // in z-Richtung
                    {
                        OrtBx[MapArray(i2), MapArray(j2), MapArray(k2)] = i2 * permanentMagnet.Bsw;
                        sx = OrtBx[MapArray(i2), MapArray(j2), MapArray(k2)];
                        OrtBy[MapArray(i2), MapArray(j2), MapArray(k2)] = j2 * permanentMagnet.Bsw;
                        sy = OrtBy[MapArray(i2), MapArray(j2), MapArray(k2)];
                        OrtBz[MapArray(i2), MapArray(j2), MapArray(k2)] = k2 * permanentMagnet.Bsw;
                        sz = OrtBz[MapArray(i2), MapArray(j2), MapArray(k2)];

                        this.BerechneHges(permanentMagnet, out Hx, out Hy, out Hz, sx, sy, sz);

                        Bx[MapArray(i2), MapArray(j2), MapArray(k2)] = Physics.Muo * Hx; // Telsa
                        By[MapArray(i2), MapArray(j2), MapArray(k2)] = Physics.Muo * Hy; // Telsa
                        Bz[MapArray(i2), MapArray(j2), MapArray(k2)] = Physics.Muo * Hz; // Telsa
                    }
                }

                Print(".");
            }

            PrintLine(" Obere Spule ist durchgerechnet.");
            //  Dann die untere Spule dazu addieren: 
            permanentMagnet.MEyo = -permanentMagnet.MEyo;   // Position der unteren Spule

            for (int i2 = -Bn; i2 < Bn; i2++)  // in x-Richtung
            {
                for (int j2 = -Bn; j2 < Bn; j2++)  // in y-Richtung
                {
                    for (int k2 = -Bn; k2 < Bn; k2++)  // in z-Richtung
                    {
                        OrtBx[MapArray(i2), MapArray(j2), MapArray(k2)] = i2 * permanentMagnet.Bsw;
                        sx = OrtBx[MapArray(i2), MapArray(j2), MapArray(k2)];
                        OrtBy[MapArray(i2), MapArray(j2), MapArray(k2)] = j2 * permanentMagnet.Bsw;
                        sy = OrtBy[MapArray(i2), MapArray(j2), MapArray(k2)];
                        OrtBz[MapArray(i2), MapArray(j2), MapArray(k2)] = k2 * permanentMagnet.Bsw;
                        sz = OrtBz[MapArray(i2), MapArray(j2), MapArray(k2)];

                        this.BerechneHges(permanentMagnet, out Hx, out Hy, out Hz, sx, sy, sz);

                        Bx[MapArray(i2), MapArray(j2), MapArray(k2)] = Bx[MapArray(i2), MapArray(j2), MapArray(k2)] + Physics.Muo * Hx; // Telsa
                        By[MapArray(i2), MapArray(j2), MapArray(k2)] = By[MapArray(i2), MapArray(j2), MapArray(k2)] + Physics.Muo * Hy; // Telsa
                        Bz[MapArray(i2), MapArray(j2), MapArray(k2)] = Bz[MapArray(i2), MapArray(j2), MapArray(k2)] + Physics.Muo * Hz; // Telsa
                    }
                }

                Print(".");
            }

            PrintLine(" Untere Spule ist durchgerechnet.");
            PrintLine();
            permanentMagnet.MEyo = -permanentMagnet.MEyo;   // MEyo zurücksetzen.
            PrintLine("Gesamtes Feld am Koordinaten-Ursprung: ");
            PrintLine(Bx[0, 0, 0] + ", " + By[0, 0, 0] + ", " + Bz[0, 0, 0] + " T");
            PrintLine();
            PrintLine("Gesamtes Feld im Zentrum der oberen Spule:");
            // oberen Spulenmittelpunkt suchen: 
            sx = 0; sy = permanentMagnet.MEyo; sz = 0;
            this.BerechneHges(permanentMagnet, out Hx, out Hy, out Hz, sx, sy, sz);
            BXmax = Physics.Muo * Hx;
            BYmax = Physics.Muo * Hy;
            BZmax = Physics.Muo * Hz;
            // unteren Spulenmittelpunkt suchen: 
            sx = 0; sy = -permanentMagnet.MEyo; sz = 0;
            this.BerechneHges(permanentMagnet, out Hx, out Hy, out Hz, sx, sy, sz);
            BXmax = BXmax + Physics.Muo * Hx;
            BYmax = BYmax + Physics.Muo * Hy;
            BZmax = BZmax + Physics.Muo * Hz;

            PrintLine(BXmax + ", " + BYmax + ", " + BZmax + " T");
            PrintLine("Ist dieses Feld gewünscht ? ? ! ? ? ! ? ?");
        }

/*
        private void Magnetfeld_anzeigen()
        {
            PrintLine("Feld 'Magnetische Induktion' des Dauermagneten:");

            for (int i = -Bn; i < Bn; i++)  // in x-Richtung
            {
                for (int j = -Bn; j < Bn; j++)  // in y-Richtung
                {
                    for (int k = -Bn; k < Bn; k++) // in z-Richtung
                    {
                        Print("x,y,z=" + OrtBx[i, j, k] * 100 + ", " + OrtBy[i, j, k] * 100 + ", " + OrtBz[i, j, k] * 100 + "cm => B=");
                        Print(Bx[MapArray(i), MapArray(j), MapArray(k)] + ", ");
                        Print(By[MapArray(i), MapArray(j), MapArray(k)] + ", ");
                        Print(Bz[MapArray(i), MapArray(j), MapArray(k)] + " T ");
                    }
                }
            }
        }
*/

        private void StromverteilungZuweisen03(PermanentMagnetParameters permanentMagnet)
        {
            PrintLine("Kontrolle der Magnetfeld-Emulations-Spulen:");
            for (int i = 0; i < MESEanz / 2; i++)
            {
                MESEx[i] = permanentMagnet.MEro * Math.Cos((i - 1) / (MESEanz / 2) * 2.0 * Math.PI);    // Orte der oberen Magnetfeld-Emulations-Spulen-Elemente
                MESEy[i] = permanentMagnet.MEyo;                                                        // Orte der oberen Magnetfeld-Emulations-Spulen-Elemente
                MESEz[i] = permanentMagnet.MEro * Math.Sin((i - 1) / (MESEanz / 2) * 2.0 * Math.PI);    // Orte der oberen Magnetfeld-Emulations-Spulen-Elemente
                MESEdx[i] = -Math.Sin((i - 1) / (MESEanz / 2) * 2.0 * Math.PI);         // Laufrichtungen der Magnetfeld-Emulations-Spulen-Elemente
                MESEdy[i] = 0;                                                          // Laufrichtungen der Magnetfeld-Emulations-Spulen-Elemente
                MESEdz[i] = Math.Cos((i - 1) / (MESEanz / 2) * 2.0 * Math.PI);          // Laufrichtungen der Magnetfeld-Emulations-Spulen-Elemente
            }

            for (int i = (MESEanz / 2) + 1; i < MESEanz; i++)
            {
                MESEx[i] = permanentMagnet.MEro * Math.Cos((i - 1) / (MESEanz / 2) * 2.0 * Math.PI);    // Orte der unteren Magnetfeld-Emulations-Spulen-Elemente
                MESEy[i] = -permanentMagnet.MEyo;                                                       // Orte der unteren Magnetfeld-Emulations-Spulen-Elemente
                MESEz[i] = permanentMagnet.MEro * Math.Sin((i - 1) / (MESEanz / 2) * 2.0 * Math.PI);    // Orte der unteren Magnetfeld-Emulations-Spulen-Elemente
                MESEdx[i] = -Math.Sin((i - 1) / (MESEanz / 2) * 2.0 * Math.PI);         // Laufrichtungen der Magnetfeld-Emulations-Spulen-Elemente
                MESEdy[i] = 0;                                                          // Laufrichtungen der Magnetfeld-Emulations-Spulen-Elemente
                MESEdz[i] = Math.Cos((i - 1) / (MESEanz / 2) * 2.0 * Math.PI);          // Laufrichtungen der Magnetfeld-Emulations-Spulen-Elemente
            }
        }

        private void SpulenZuweisen() // Spule für den Input der Steuer-Energie
        {
            // Zuerst die Stützpunkte des Polygone:
            for (int i = 0; i < 2 * zo; i++)
            {   // Anfang äußerst links unten, gehe erst in z-Richtung
                SpIx[i + 1] = -xo * Spsw;
                SpIy[i + 1] = -yo * Spsw;
                SpIz[i + 1] = (i - zo) * Spsw;     // Stützpunkt

                SpTx[i + 1] = +xo * Spsw;
                SpTy[i + 1] = -yo * Spsw;
                SpTz[i + 1] = (i - zo) * Spsw;     // Stützpunkt

                SIx[i + 1] = -xo * Spsw;
                SIy[i + 1] = -yo * Spsw;
                SIz[i + 1] = (0.5 + i - zo) * Spsw; // Ort,Mittelpunkt

                STx[i + 1] = -xo * Spsw;
                STy[i + 1] = -yo * Spsw;
                STz[i + 1] = (0.5 + i - zo) * Spsw; // Ort,Mittelpunkt

                dSIx[i + 1] = 0;
                dSIy[i + 1] = 0;
                dSIz[i + 1] = +Spsw;           // Richtungsvektor

                dSTx[i + 1] = 0;
                dSTy[i + 1] = 0;
                dSTz[i + 1] = +Spsw;           // Richtungsvektor
            }

            for (int i = 0; i < 2 * yo; i++)
            {   // gehe dann weiter in y-Richtung
                SpIx[2 * zo + i + 1] = -xo * Spsw;
                SpIy[2 * zo + i + 1] = (i - yo) * Spsw;
                SpIz[2 * zo + i + 1] = +zo * Spsw;  // Stützpunkt

                SpTx[2 * zo + i + 1] = +xo * Spsw;
                SpTy[2 * zo + i + 1] = (i - yo) * Spsw;
                SpTz[2 * zo + i + 1] = +zo * Spsw;  // Stützpunkt

                SIx[2 * zo + i + 1] = -xo * Spsw;
                SIy[2 * zo + i + 1] = (0.5 + i - yo) * Spsw;
                SIz[2 * zo + i + 1] = +zo * Spsw;  // Ort,Mittelpunkt

                STx[2 * zo + i + 1] = +xo * Spsw;
                STy[2 * zo + i + 1] = (0.5 + i - yo) * Spsw;
                STz[2 * zo + i + 1] = +zo * Spsw;  // Ort,Mittelpunkt

                dSIx[2 * zo + i + 1] = 0;
                dSIy[2 * zo + i + 1] = Spsw;
                dSIz[2 * zo + i + 1] = 0;         // Richtungsvektor

                dSTx[2 * zo + i + 1] = 0;
                dSTy[2 * zo + i + 1] = Spsw;
                dSTz[2 * zo + i + 1] = 0;         // Richtungsvektor
            }

            for (int i = 0; i < 2 * zo; i++)
            {   // gehe dann in z-Richtung zurück
                SpIx[2 * zo + 2 * yo + i + 1] = -xo * Spsw;
                SpIy[2 * zo + 2 * yo + i + 1] = yo * Spsw;
                SpIz[2 * zo + 2 * yo + i + 1] = (zo - i) * Spsw;     // Stützpunkt

                SpTx[2 * zo + 2 * yo + i + 1] = +xo * Spsw;
                SpTy[2 * zo + 2 * yo + i + 1] = yo * Spsw;
                SpTz[2 * zo + 2 * yo + i + 1] = (zo - i) * Spsw;     // Stützpunkt

                SIx[2 * zo + 2 * yo + i + 1] = -xo * Spsw;
                SIy[2 * zo + 2 * yo + i + 1] = yo * Spsw;
                SIz[2 * zo + 2 * yo + i + 1] = (zo - i - 0.5) * Spsw; // Ort,Mittelpunkt

                STx[2 * zo + 2 * yo + i + 1] = +xo * Spsw;
                STy[2 * zo + 2 * yo + i + 1] = yo * Spsw;
                STz[2 * zo + 2 * yo + i + 1] = (zo - i - 0.5) * Spsw; // Ort,Mittelpunkt

                dSIx[2 * zo + 2 * yo + i + 1] = 0;
                dSIy[2 * zo + 2 * yo + i + 1] = 0;
                dSIz[2 * zo + 2 * yo + i + 1] = -Spsw;           // Richtungsvektor

                dSTx[2 * zo + 2 * yo + i + 1] = 0;
                dSTy[2 * zo + 2 * yo + i + 1] = 0;
                dSTz[2 * zo + 2 * yo + i + 1] = -Spsw;           // Richtungsvektor
            }

            for (int i = 0; i < 2 * yo; i++)
            {   // und zu guter Letzt in y-Richtung wieder runter
                SpIx[4 * zo + 2 * yo + i + 1] = -xo * Spsw;
                SpIy[4 * zo + 2 * yo + i + 1] = (yo - i) * Spsw;
                SpIz[4 * zo + 2 * yo + i + 1] = -zo * Spsw; // Stützpunkt

                SpTx[4 * zo + 2 * yo + i + 1] = +xo * Spsw;
                SpTy[4 * zo + 2 * yo + i + 1] = (yo - i) * Spsw;
                SpTz[4 * zo + 2 * yo + i + 1] = -zo * Spsw; // Stützpunkt

                SIx[4 * zo + 2 * yo + i + 1] = -xo * Spsw;
                SIy[4 * zo + 2 * yo + i + 1] = (yo - i - 0.5) * Spsw;
                SIz[4 * zo + 2 * yo + i + 1] = -zo * Spsw; // Ort,Mittelpunkt

                STx[4 * zo + 2 * yo + i + 1] = +xo * Spsw;
                STy[4 * zo + 2 * yo + i + 1] = (yo - i - 0.5) * Spsw;
                STz[4 * zo + 2 * yo + i + 1] = -zo * Spsw; // Ort,Mittelpunkt

                dSIx[4 * zo + 2 * yo + i + 1] = 0;
                dSIy[4 * zo + 2 * yo + i + 1] = -Spsw;
                dSIz[4 * zo + 2 * yo + i + 1] = 0;        // Richtungsvektor

                dSTx[4 * zo + 2 * yo + i + 1] = 0;
                dSTy[4 * zo + 2 * yo + i + 1] = -Spsw;
                dSTz[4 * zo + 2 * yo + i + 1] = 0;        // Richtungsvektor
            }   // Der letzte Punkt ist dem ersten identisch

            SpN = 4 * zo + 4 * yo + 1;

            PrintLine("Anzahl der Punkte der Spulen-Linienaufteilung: von 1 - " + SpN);
            if (SpN > SpNmax)
            {
                PrintLine("--- ERROR --- zu viele Spulen-Linienelemente");
                PrintLine("--- ABHILFE -> Array grösser dimensionieren");
                Environment.Exit(0);
            }

            // Dann die Flächenelemente:
            for (int j = 1; j < 2 * yo; j++)
            {
                for (int i = 1; i < 2 * zo; i++)
                {
                    FlIx[i + (j - 1) * 2 * zo] = -xo * Spsw;
                    FlIy[i + (j - 1) * 2 * zo] = (j - 0.5 - yo) * Spsw;
                    FlIz[i + (j - 1) * 2 * zo] = (i - 0.5 - zo) * Spsw;

                    FlTx[i + (j - 1) * 2 * zo] = +xo * Spsw;
                    FlTy[i + (j - 1) * 2 * zo] = (j - 0.5 - yo) * Spsw;
                    FlTz[i + (j - 1) * 2 * zo] = (i - 0.5 - zo) * Spsw;
                }
            }

            FlN = 4 * zo * yo;

            PrintLine("Anzahl der Flächenelemente jeder Spulen: von 1 - " + FlN);
            if (FlN > FlNmax)
            {
                PrintLine("--- ERROR --- zu viele Spulen-Flächenelemente");
                PrintLine("--- ABHILFE -> Array grösser dimensionieren");
                Environment.Exit(0);
            }
        }

/*
        private void Spulen_anzeigen()  // Spule für den Input der Steuer-Energie
        {
            PrintLine("Input-Sp.-> Stützpunkte des Polygons, Orte der FE, Richtungsvektoren der FE:");
            for (int i = 1; i < SpN; i++)
            {
                PrintLine("SP [" + i + "]= " + SpIx[i] * 100 + ", " + SpIy[i] * 100 + ", " + SpIz[i] * 100 + " cm ");
                PrintLine("ORT[" + i + "]= " + SIx[i] * 100 + ", " + SIy[i] * 100 + ", " + SIz[i] * 100 + " cm ");
                PrintLine("RV [" + i + "]= " + dSIx[i] * 100 + ", " + dSIy[i] * 100 + ", " + dSIz[i] * 100 + " cm ");
            }

            PrintLine("Turbo-Sp.-> Stützpunkte des Polygons, Orte der FE, Richtungsvektoren der FE:");
            for (int i = 1; i < SpN; i++)
            {
                PrintLine("SP [" + i + "]= " + SpTx[i] * 100 + ", " + SpTy[i] * 100 + ", " + SpTz[i] * 100 + " cm ");
                PrintLine("ORT[" + i + "]= " + STx[i] * 100 + ", " + STy[i] * 100 + ", " + STz[i] * 100 + " cm ");
                PrintLine("RV [" + i + "]= " + dSTx[i] * 100 + ", " + dSTy[i] * 100 + ", " + dSTz[i] * 100 + " cm ");
            }

            PrintLine("Input-Spule -> Flächenelemente, deren Mittelpunktspositionen:");
            for (int i = 1; i < FlN; i++)
            {
                Print("x,y,z[" + i + "]= " + FlIx[i] * 100 + ", " + FlIy[i] * 100 + ", " + FlIz[i] * 100 + " cm ");
            }

            PrintLine("Turbo-Spule -> Flächenelemente, deren Mittelpunktspositionen:");
            for (int i = 1; i < FlN; i++)
            {
                Print("x,y,z[" + i + "]= " + FlTx[i] * 100 + ", " + FlTy[i] * 100 + ", " + FlTz[i] * 100 + " cm ");
            }

            PrintLine("--------------------------------------------------------------------------------");
        }
*/

        private void MagnetDrehen(double fi) // Drehen um Drehwinkel "fi":
        {
            fi = fi / 180 * Math.PI;  // Umrechnen des Winkels in Radianten
            for (int i = -Bn; i < Bn; i++)  // x-Anteile
            {
                for (int j = -Bn; j < Bn; j++)  // y-Anteile
                {
                    for (int k = -Bn; k < Bn; k++)  // z-Anteile
                    {
                        // Drehung der Ortsvektoren:
                        OrtBxDR[MapArray(i), MapArray(j), MapArray(k)] = +OrtBx[MapArray(i), MapArray(j), MapArray(k)] * Math.Cos(-fi) + OrtBy[MapArray(i), MapArray(j), MapArray(k)] * Math.Sin(-fi);
                        OrtByDR[MapArray(i), MapArray(j), MapArray(k)] = -OrtBx[MapArray(i), MapArray(j), MapArray(k)] * Math.Sin(-fi) + OrtBy[MapArray(i), MapArray(j), MapArray(k)] * Math.Cos(-fi);
                        OrtBzDR[MapArray(i), MapArray(j), MapArray(k)] = +OrtBz[MapArray(i), MapArray(j), MapArray(k)];
                        // Drehung der Feldstärke-Vektoren:
                        BxDR[MapArray(i), MapArray(j), MapArray(k)] = +Bx[MapArray(i), MapArray(j), MapArray(k)] * Math.Cos(-fi) + By[MapArray(i), MapArray(j), MapArray(k)] * Math.Sin(-fi);
                    }
                }
            }
        }

        private void FeldstaerkeAmOrtSuchen(double xpos, double ypos, double zpos) // an dem Ort suche ich die Feldstärke
        {
            // Zuerst suche ich, welcher Feldort zum xpos,ypos,zpos den kürzesten Abstand hat.
            int ixo = 0; int iyo = 0; int izo = 0;

            double disto = Math.Sqrt(Physics.Sqr(xpos - this.OrtBxDR[ixo, iyo, izo]) + Physics.Sqr(ypos - this.OrtByDR[ixo, iyo, izo]) + Physics.Sqr(zpos - this.OrtBzDR[ixo, iyo, izo]));

            //  PrintLine("Anfangs-Abstand vom Nullpunkt: ",disto*100:1:15," cm"); 
            for (int ix = -Bn; ix < Bn; ix++)  // x-Suche
            {
                for (int iy = -Bn; iy < Bn; iy++)  // y-Suche
                {
                    for (int iz = -Bn; iz < Bn; iz++)  // z-Suche
                    {
                        double dist = Math.Sqrt(Physics.Sqr(xpos - OrtBxDR[MapArray(ix), MapArray(iy), MapArray(iz)]) + Physics.Sqr(ypos - OrtByDR[MapArray(ix), MapArray(iy), MapArray(iz)]) + Physics.Sqr(zpos - OrtBzDR[MapArray(ix), MapArray(iy), MapArray(iz)]));

                        if (dist <= disto)
                        {
                            ixo = ix; iyo = iy; izo = iz;
                            disto = dist;
                        }
                    }
                }
            }

            //  PrintLine("Punkt Nummer (ixo,iyo,izo): ",ixo,",  ",iyo,",  ",izo); 
            // Dann gebe ich das Magnetfeld eben)rt an:
            //  PrintLine("Magnetfeld )rt: ",BxDR[ixo,iyo,izo]:8:4,", ",ByDR[ixo,iyo,izo]:8:4,", ",BzDR[ixo,iyo,izo]:8:4," T "); 
            // Jetzt brauche ich noch den magnetischen Fluß durch das )rtige Spulen-Flächenelement:
            PsiSFE = BxDR[MapArray(ixo), MapArray(iyo), MapArray(izo)] * Spsw * Spsw;      // nach *1 von S.3
            //  PrintLine("Magn. Fluss durch Spulen-Flächenelement: ",PsiSFE," T*m^2"); 
        }

        private void GesamtflussDurchInputSpule()  // gemäß *2 von S.3
        {
            PsiGES = 0;
            for (int i = 1; i < FlN; i++)
            {
                this.FeldstaerkeAmOrtSuchen(FlIx[i], FlIy[i], FlIz[i]);
                PsiGES = PsiGES + PsiSFE;
            }
        }

        private void GesamtflussDurchTurboSpule()  // gemäß *2 von S.3
        {
            PsiGES = 0;
            for (int i = 1; i < FlN; i++)
            {
                this.FeldstaerkeAmOrtSuchen(FlTx[i], FlTy[i], FlTz[i]);
                PsiGES = PsiGES + PsiSFE;
            }
        }

        private double FlussI(double alpha)
        {
            double merk = B1I * Math.Sin(1.0 * alpha);
            merk = merk + B2I * Math.Sin(2.0 * alpha);
            merk = merk + B3I * Math.Sin(3.0 * alpha);
            merk = merk + B4I * Math.Sin(4.0 * alpha);
            merk = merk + B5I * Math.Sin(5.0 * alpha);
            return merk;
        }

        private double FlussT(double alpha)
        {
            double merk = B1T * Math.Sin(1.0 * alpha);// Hier ist alpha in "radianten" anzugeben.
            merk = merk + B2T * Math.Sin(2.0 * alpha);
            merk = merk + B3T * Math.Sin(3.0 * alpha);
            merk = merk + B4T * Math.Sin(4.0 * alpha);
            merk = merk + B5T * Math.Sin(5.0 * alpha);
            return merk;
        }

        private double SchnellDrehmoment(double winkel)
        {
            return B1dreh * Math.Sin(winkel + phase);
        }

        private void MagfeldTurboBerechnen(SimulationParameters param, double rx, double ry, double rz, double Strom)
        {
            double sx, sy, sz;    // Orte der Leiterschleifen-Elemente
            double dsx, dsy, dsz;    // Laufrichtungen der Leiterschleifen-Elemente
            int AnzLSE;   // Anzahl der Leiterschleifen-Elemente
            double smrx, smry, smrz; // Differenzen für das nachfolgende Kreuzprodukt
            double krpx, krpy, krpz; // Kreuzprodukt in Biot-Savart
            double smrbetrhoch3; // Betragsbildung für Nenner
            double dHx, dHy, dHz; // Infinitesimales Magnetfeld
            double Hgesx, Hgesy, Hgesz;// Gesamt-Magnetfeld der Input-Spule am Aufpunkt

            AnzLSE = SpN - 1;
            if (AnzLSE != 4 * yo + 4 * zo)
            {
                PrintLine("Da timmt wat nich: Vernetzung der felderzeugenden Spule ist falsch.");
                PrintLine("Problem bei der TURBO-Spule:");
                PrintLine("Anzahl der Stützpunkte der Spulen, AnzLSE = " + AnzLSE);
                PrintLine("Hingegen: 4*yo+4*zo = " + 4 * yo + 4 * zo);
                Environment.Exit(1);
            }

            // Ort und Richtungsvektoren der Leiterschleifen-Elemente feststellen, Feld am Aufpunkt ausrechnen nach Biot-Savart:
            Hgesx = 0; Hgesy = 0; Hgesz = 0;
            for (int i = 1; i < AnzLSE; i++)
            {
                sx = SpTx[i]; sy = SpTy[i]; sz = SpTz[i];     // Orte der Leiterschleifen-Elemente
                dsx = dSTx[i]; dsy = dSTy[i]; dsz = dSTz[i];    // Richtungsvektoren der Leiterschleifen-Elemente
                smrx = sx - rx; smry = sy - ry; smrz = sz - rz;   // Differenzen für das nachfolgende Kreuzprodukt
                krpx = dsy * smrz - dsz * smry; krpy = dsz * smrx - dsx * smrz; krpz = dsx * smry - dsy * smrx; // Kreuzprodukt
                smrbetrhoch3 = Math.Sqrt(Physics.Sqr(smrx) + Physics.Sqr(smry) + Physics.Sqr(smrz));
                if (smrbetrhoch3 < Spsw / 1000)
                {
                    PrintLine("Mechanische Kollision -> Magnet berührt Turbo-Spule. STOP.");
                    PrintLine("Spulen-Element bei: " + sx + ", " + sy + ", " + sz + "m.");
                    PrintLine("Magnet-Ort     bei: " + rx + ", " + ry + ", " + rz + "m.");
                    Environment.Exit(1);
                }

                smrbetrhoch3 = smrbetrhoch3 * smrbetrhoch3 * smrbetrhoch3;  // Betragsbildung für Nenner in Biot-Savart
                dHx = Strom * krpx / 4 / Math.PI / smrbetrhoch3;           // Finites Magnetfeld des Leiterschleifen-Elements
                dHy = Strom * krpy / 4 / Math.PI / smrbetrhoch3;
                dHz = Strom * krpz / 4 / Math.PI / smrbetrhoch3;
                Hgesx = Hgesx + dHx; Hgesy = Hgesy + dHy; Hgesz = Hgesz + dHz; // Summation der Feldelemente
            }       // Vorzeichen in der nachfolgenden Zeile gemäß technischer Stromrichtung.

            BTx = -Physics.Muo * Hgesx * param.Inductor.Nturbo;
            BTy = -Physics.Muo * Hgesy * param.Inductor.Nturbo;
            BTz = -Physics.Muo * Hgesz * param.Inductor.Nturbo;
        }

        private void MagfeldInputBerechnen(SimulationParameters param, double rx, double ry, double rz, double Strom)
        {
            double sx, sy, sz;    // Orte der Leiterschleifen-Elemente
            double dsx, dsy, dsz;    // Laufrichtungen der Leiterschleifen-Elemente
            int AnzLSE;   // Anzahl der Leiterschleifen-Elemente
            double smrx, smry, smrz; // Differenzen für das nachfolgende Kreuzprodukt
            double krpx, krpy, krpz; // Kreuzprodukt in Biot-Savart
            double smrbetrhoch3; // Betragsbildung für Nenner
            double dHx, dHy, dHz; // Infinitesimales Magnetfeld
            double Hgesx, Hgesy, Hgesz;// Gesamt-Magnetfeld der Input-Spule am Aufpunkt

            AnzLSE = SpN - 1;
            if (AnzLSE != 4 * yo + 4 * zo)
            {
                PrintLine("Da timmt wat nich: Vernetzung der felderzeugenden Spule ist falsch.");
                PrintLine("Problem bei der INPUT-Spule:");
                PrintLine("Anzahl der Stützpunkte der Spulen, AnzLSE = " + AnzLSE);
                PrintLine("Hingegen: 4*yo+4*zo = " + 4 * yo + 4 * zo);
                Environment.Exit(1);
            }

            // Ort und Richtungsvektoren der Leiterschleifen-Elemente feststellen, Feld am Aufpunkt ausrechnen nach Biot-Savart:
            Hgesx = 0; Hgesy = 0; Hgesz = 0;
            for (int i = 1; i < AnzLSE; i++)
            {
                sx = SpIx[i]; sy = SpIy[i]; sz = SpIz[i];     // Orte der Leiterschleifen-Elemente
                dsx = dSIx[i]; dsy = dSIy[i]; dsz = dSIz[i];    // Richtungsvektoren der Leiterschleifen-Elemente
                smrx = sx - rx; smry = sy - ry; smrz = sz - rz;   // Differenzen für das nachfolgende Kreuzprodukt
                krpx = dsy * smrz - dsz * smry; krpy = dsz * smrx - dsx * smrz; krpz = dsx * smry - dsy * smrx; // Kreuzprodukt
                smrbetrhoch3 = Math.Sqrt(Physics.Sqr(smrx) + Physics.Sqr(smry) + Physics.Sqr(smrz));
                if (smrbetrhoch3 < Spsw / 1000)
                {
                    PrintLine("Mechanische Kollision -> Magnet berührt Input-Spule. STOP.");
                    PrintLine("Spulen-Element bei: " + sx + ", " + sy + ", " + sz + "m.");
                    PrintLine("Magnet-Ort     bei: " + rx + ", " + ry + ", " + rz + "m.");
                    Environment.Exit(1);
                }
                smrbetrhoch3 = smrbetrhoch3 * smrbetrhoch3 * smrbetrhoch3;  // Betragsbildung für Nenner in Biot-Savart
                dHx = Strom * krpx / 4.0 / Math.PI / smrbetrhoch3;           // Finites Magnetfeld des Leiterschleifen-Elements
                dHy = Strom * krpy / 4.0 / Math.PI / smrbetrhoch3;
                dHz = Strom * krpz / 4.0 / Math.PI / smrbetrhoch3;
                Hgesx = Hgesx + dHx; Hgesy = Hgesy + dHy; Hgesz = Hgesz + dHz; // Summation der Feldelemente
            }       // Vorzeichen in der nachfolgenden Zeile gemäß technischer Stromrichtung.

            BIx = -Physics.Muo * Hgesx * param.Inductor.Ninput; BIy = -Physics.Muo * Hgesy * param.Inductor.Ninput; BIz = -Physics.Muo * Hgesz;
        }

        private double Drehmoment(SimulationParameters param, double alpha) // Das Argument ist der Winkel der Magnetstellung "alpha"
        {
            double Idlx, Idly, Idlz; // Kartesische Komponenten von dl-Vektor nach (*1 von S.11)
            double Bxlok, Bylok, Bzlok; // lokale Magnetfeld-Werte
            double FLx, FLy, FLz; // Lorentz-Kraft als Kreuzprodukt
            double dMx, dMy, dMz; // Drehmoment, das jedes Leiterschleifen-Element auf seinen gesamten Magneten ausübt.
            double MgesX, MgesY, MgesZ; // Gesamt-Drehmoment der Turbo- & Input- Spule auf den Magneten (aus Emulations-Spulen).
            double rx, ry, rz; // Ortsangabe der Magnetfeld-Emulationsspulen-Elemente nach Dreh-Transformation

            MgesX = 0; MgesY = 0; MgesZ = 0;
            for (int i = 1; i < MESEanz; i++)
            {
                // Wir {nen mit der Berechnung der Lorentz-Kraft auf jedes einzelne Element der Magnetfeld-Emulations-Spule
                Idlx = param.PermanentMagnet.MEI * MESEdx[i] * 4.0 * Math.PI * param.PermanentMagnet.MEro / MESEanz; // Magnetfeld-Emulations-Leiterschleifen-Element
                Idly = param.PermanentMagnet.MEI * MESEdy[i] * 4.0 * Math.PI * param.PermanentMagnet.MEro / MESEanz; // Magnetfeld-Emulations-Leiterschleifen-Element
                Idlz = param.PermanentMagnet.MEI * MESEdz[i] * 4.0 * Math.PI * param.PermanentMagnet.MEro / MESEanz; // Magnetfeld-Emulations-Leiterschleifen-Element
                // Es folgt die Berechnung der Magnetfeld-Stärke der Input- und Turbo- Spule am Ort des Leiterschleifen-Elements
                this.MagfeldInputBerechnen(param, MESEx[i], MESEy[i], MESEz[i], qpoI); // Strom durch die Input-Spule einstellen !!
                this.MagfeldTurboBerechnen(param, MESEx[i], MESEy[i], MESEz[i], qpoT); // Strom "QP" durch die Input-Spule einstellen !!
                Bxlok = BIx + BTx; Bylok = BIy + BTy; Bzlok = BIz + BTz;   // lokales B-feld der Input- und der Turbo- Spule am Ort des Leiterschleifen-Elements
                // Kreuzprodukt bilden zur Berechnung der Lorentz-Kraft:
                FLx = Idly * Bzlok - Idlz * Bylok;
                FLy = Idlz * Bxlok - Idlx * Bzlok;
                FLz = Idlx * Bylok - Idly * Bxlok;
                // Kontrolle der Lorentz-Kraft:
                /*
                    PrintLine("Ort: " + MESEx[i] + ", " + MESEy[i] + ", " + MESEz[i]);
                    PrintLine(" dl: " + MESEdx[i] + ", " + MESEdy[i] + ", " + MESEdz[i]);
                    PrintLine("FLo: " + FLx + ", " + FLy + ", " + FLz);
                */
                // Den wirkenden Ort, an dem das Drehmoment angreift, stellen wir durch eine Drehtransformation gemäß *1 von S.12 fest:
                rx = +MESEx[i] * Math.Cos(-alpha) + MESEy[i] * Math.Sin(-alpha);
                ry = -MESEx[i] * Math.Sin(-alpha) + MESEy[i] * Math.Cos(-alpha);
                rz = MESEz[i];
                // Daraus berechnen wir nun das zugehörige Drehmoment-Element, das dieses Lorenzt-Kraft-Element auf dem Magneten ausübt:
                dMx = ry * FLz - rz * FLy;           // Drehmoment als Kreuzprodukt M = r x F 
                dMy = rz * FLx - rx * FLz;
                dMz = rx * FLy - ry * FLx;
                // Kontrolle des Drehmoments:
                //    PrintLine("Dreh:",dMx,", ",dMy,", ",dMz); Wait;         
                MgesX = MgesX + dMx; // Summation aller einzelnen Drehmoment-Elemente zum Gesamt-Drehmoment.
                MgesY = MgesY + dMy; // in drei kartesischen Komponenten
                MgesZ = MgesZ + dMz; // Wegen der Lagerung des Magneten sMath.PIelt nur die z-Komponente des Drehmoments mit.
            }                // Der Magnet hat nämlich eine starre Achse und rotiert nur um die z-Achse.
            //  PrintLine("Drehmoment:",MgesX:20,", ",Mgesy:20,", ",Mgesz:20); 
            return MgesZ;
        }

        private void DatenSpeichern(SimulationParameters param)
        {
            using (var stream = File.OpenWrite("schonda"))
            using (var fout = new StreamWriter(stream))
            {
                // Zuerst die Parameter:
                fout.WriteLine(Spsw);
                fout.WriteLine(xo);
                fout.WriteLine(yo);
                fout.WriteLine(zo);
                fout.WriteLine(param.Inductor.Ninput);
                fout.WriteLine(param.Inductor.Nturbo);
                fout.WriteLine(param.PermanentMagnet.Bsw);
                fout.WriteLine(param.PermanentMagnet.MEyo);
                fout.WriteLine(param.PermanentMagnet.MEro);
                fout.WriteLine(param.PermanentMagnet.MEI);

                // Dann das Magnetfeld:}   {Die Schritt-Anzahl Bn bleibt "Const."
                for (int i = -Bn; i < Bn; i++)  // in x-Richtung
                {
                    for (int j = -Bn; j < Bn; j++)  // in y-Richtung
                    {
                        for (int k = -Bn; k < Bn; k++)  // in z-Richtung
                        {
                            fout.WriteLine(OrtBx[MapArray(i), MapArray(j), MapArray(k)]);
                            fout.WriteLine(OrtBy[MapArray(i), MapArray(j), MapArray(k)]);
                            fout.WriteLine(OrtBz[MapArray(i), MapArray(j), MapArray(k)]);
                            fout.WriteLine(Bx[MapArray(i), MapArray(j), MapArray(k)]);
                            fout.WriteLine(By[MapArray(i), MapArray(j), MapArray(k)]);
                            fout.WriteLine(Bz[MapArray(i), MapArray(j), MapArray(k)]);
                        }
                    }
                }

                // Spulenzuweisung und Stromverteilung brauche ich nicht speichern, die kann man rechnen lassen.
                // Die Drehmoments-Parameter muß ich auch abspeichern:
                fout.WriteLine(B1T);
                fout.WriteLine(B2T);
                fout.WriteLine(B3T);
                fout.WriteLine(B4T);
                fout.WriteLine(B5T);
                fout.WriteLine(B1I);
                fout.WriteLine(B2I);
                fout.WriteLine(B3I);
                fout.WriteLine(B4I);
                fout.WriteLine(B5I);
                fout.WriteLine(B1dreh);
                fout.WriteLine(phase);
                fout.WriteLine("Die Daten sind alle herausgeschrieben.");
            }
        }

        private bool AlteParameterVergleichen(SimulationParameters param)
        {
            double x;  // Parameter zum Einlesen
            int n; // Parameter zum Einlesen
            bool schonda = true;

            if (!File.Exists("schonda"))
            {
                return false;
            }

            using (var stream = File.OpenRead("schonda"))
            using (var fin = new StreamReader(stream))
            {
                // Zuerst die Parameter:
                x = double.Parse(fin.ReadLine());
                if (x != Spsw) { schonda = false; PrintLine(" Spsw  geändert"); }
                n = int.Parse(fin.ReadLine());
                if (n != xo) { schonda = false; PrintLine("  xo   geändert"); }
                n = int.Parse(fin.ReadLine());
                if (n != yo) { schonda = false; PrintLine("  yo   geändert"); }
                n = int.Parse(fin.ReadLine());
                if (n != zo) { schonda = false; PrintLine("  zo   geändert"); }
                n = int.Parse(fin.ReadLine());
                if (n != param.Inductor.Ninput) { schonda = false; PrintLine("Ninput geändert"); }
                n = int.Parse(fin.ReadLine());
                if (n != param.Inductor.Nturbo) { schonda = false; PrintLine("Nturbo geändert"); }
                x = double.Parse(fin.ReadLine());
                if (x != param.PermanentMagnet.Bsw) { schonda = false; PrintLine("  Bsw  geändert"); }
                x = double.Parse(fin.ReadLine());
                if (x != param.PermanentMagnet.MEyo) { schonda = false; PrintLine(" MEyo  geändert"); }
                x = double.Parse(fin.ReadLine());
                if (x != param.PermanentMagnet.MEro) { schonda = false; PrintLine(" MEro  geändert"); }
                x = double.Parse(fin.ReadLine());
                if (x != param.PermanentMagnet.MEI) { schonda = false; PrintLine("  MEI  geändert"); }

                if (schonda)
                {
                    PrintLine("Die Parameter sind bereits bekannt.");
                }

                if (!schonda)
                {
                    PrintLine("Die Parameter sind neu. Es beginnt eine neue Vernetzung.");
                }

                // Dann das Magnetfeld:}   {Die Schritt-Anzahl Bn bleibt "Const."
                for (int i = -Bn; i < Bn; i++)  // in x-Richtung
                {
                    for (int j = -Bn; j < Bn; j++)  // in y-Richtung
                    {
                        for (int k = -Bn; k < Bn; k++)  // in z-Richtung
                        {
                            OrtBx[MapArray(i), MapArray(j), MapArray(k)] = double.Parse(fin.ReadLine());
                            OrtBy[MapArray(i), MapArray(j), MapArray(k)] = double.Parse(fin.ReadLine());
                            OrtBz[MapArray(i), MapArray(j), MapArray(k)] = double.Parse(fin.ReadLine());
                            Bx[MapArray(i), MapArray(j), MapArray(k)] = double.Parse(fin.ReadLine());
                            By[MapArray(i), MapArray(j), MapArray(k)] = double.Parse(fin.ReadLine());
                            Bz[MapArray(i), MapArray(j), MapArray(k)] = double.Parse(fin.ReadLine());
                        }
                    }
                }

                PrintLine("Das Magnetfeld ist gelesen.");

                // Spulenzuweisung und Stromverteilung brauche ich nicht speichern, die kann man rechnen lassen.
                // Die Drehmoments-Parameter muß ich auch abspeichern:
                B1T = double.Parse(fin.ReadLine());
                B2T = double.Parse(fin.ReadLine());
                B3T = double.Parse(fin.ReadLine());
                B4T = double.Parse(fin.ReadLine());
                B5T = double.Parse(fin.ReadLine());
                B1I = double.Parse(fin.ReadLine());
                B2I = double.Parse(fin.ReadLine());
                B3I = double.Parse(fin.ReadLine());
                B4I = double.Parse(fin.ReadLine());
                B5I = double.Parse(fin.ReadLine());
                PrintLine("Die Parameter zur Berechnung des magnetischen Flusses sind gelesen.");
                B1dreh = double.Parse(fin.ReadLine());
                phase = double.Parse(fin.ReadLine());
                PrintLine("Die Parameter zur schnellen Berechnung des Drehmoments sind gelesen.");
                PrintLine("Damit steht alles für den DFEM-Algorithmus bereit.");
            }

            return schonda;
        }

        private double U7(ulong i, double iumk)       // Input-Spannung für den Input-Schwingkreis
        {
            double UAmpl;       // Spannungs-Amplitude
            ulong Pulsdauer;  // Pulsdauer in Zeitschritten von "dt"
            double Phasenshift; // Phasendifferenz zwischen oberem Umkehrpunkt und Spannungs-Impuls
            double Umerk;       // Merk-Wert für die Ausgabe der Spannung

            Umerk = 0;          // Initialisierung des Merk-Werts für die Spannungs-Ausgabe
            UAmpl = 6E-6;       // Volt, Spannungs-Amplitude
            Pulsdauer = 20;     // Pulsdauer, Anzahl der Zeitschritte von "dt"
            Phasenshift = 000;  // Phasendifferenz zwischen oberem Umkehrpunkt und Spannungs-Impuls in Zeitschritten von "dt"
            //  if( i<=Pulsdauer ) Umerk:=UAmpl;          //falls gewünscht: Start-Impuls geben
            if (i >= Pulsdauer)        // Getriggerte Pulse im Betrieb geben
            {
                if ((i >= iumk + Phasenshift) && (i <= iumk + Pulsdauer + Phasenshift)) // Hier wird das Trigger-Signal am oberen Umkehrpunkt festgemacht.
                {
                    Umerk = UAmpl;
                } // Spannung anlegen                     //Alternativ könnte man ihn z.B. auch am Nulldurchgang festmachen.
            }

            return Umerk * 0; // Wir wollen jetzt keine Energie-Zufuhr. Die Maschine soll ein Selbstläufer werden.
        }

        private double ReibungNachregeln(GeneralParameters param)
        {  // Eine kleine Schalt-Hysterese muß ich einbauen:
            double merk = cr; // Falls ich nicht außerhalb der Schalt-Hysterese liege.
            if ((phipo / 2.0 / Math.PI * 60) > 1.000001 * param.phipZiel) merk = cr * 1.000001; // Wenn's zu schnell läuft => abbremsen
            if ((phipo / 2.0 / Math.PI * 60) < 0.999999 * param.phipZiel) merk = cr * 0.999999; // Wenn's zu langsam läuft => weniger Reibung
            if ((merk < 0.8 * param.crAnfang)) merk = 0.8 * param.crAnfang; // Regelung nicht zu arg schwingen lassen, vor Allen nicht Aufschwingen.
            if ((merk > 1.2 * param.crAnfang)) merk = 1.2 * param.crAnfang; // Regelung nicht zu arg schwingen lassen, vor Allen nicht Aufschwingen.
            return merk;
        }

        public SimulationResult Calculate(SimulationParameters param)
        {
            this.Cancel = false;

            //  Eingabe-Daten-Anmerkung: Die Input-Spannung für die Input-Spule steht als letztes Unterprogramm vor dem {n des Hauptprogramms.
            //  Initialisierung - Vorgabe der Werte: }        {Wir arbeiten in SI-Einheiten
            PrintLine("DFEM-Simulation des EMDR-Motors.");
            //  Zum Lösen der Dgl. und zur Darstellung der Ergebnisse:

            //  Die beiden Spulen, vgl. Zeichnung *2 von S.1 :} {Die Spulen werden nach Vorgabe der Geometrieparameter automatisch vernetzt
            Spsw = 0.01; // Angabe in Metern: Die Spulen-Aufgliederung ist in 0.01-Meter-Schritten
            xo = 0; yo = 6; zo = 5; // Angaben in Vielfachen von Spsw} {Geometrieparameter nach Zeichnung*2 von S.1

            this.SpulenZuweisen();     // Spule für den Input der Steuer-Energie

            if (param.Inductor.Nebeninput * param.Inductor.Ueberinput != param.Inductor.Ninput)
            {
                PrintLine();
                PrintLine("Windungszahl falsch: So kann man die Input-Spule nicht anordnen !");
                Environment.Exit(1);
            }

            if (param.Inductor.Nebenturbo * param.Inductor.Ueberturbo != param.Inductor.Nturbo)
            {
                PrintLine();
                PrintLine("Windungszahl falsch: So kann man die Turbo-Spule nicht anordnen !");
                Environment.Exit(1);
            }

            bool schonda = this.AlteParameterVergleichen(param);
            if (!schonda)
            {
                this.MagnetfeldZuweisen03(param.PermanentMagnet);  // Magnetfeld zuweisen und anzeigen
            }

            this.StromverteilungZuweisen03(param.PermanentMagnet); // Stromverteilung in den Magnetfeld-Emulations-Spulen zuweisen.

            //  Abgeleitete Parameter. Die Größen werden aus den obigen Parametern berechnet, es ist keine Eingabe möglich:
            DLI = 4.0 * (yo + zo) * Spsw * param.Inductor.Ninput; // Meter} {Länge des Spulendrahtes, Input-Spule
            DLT = 4.0 * (yo + zo) * Spsw * param.Inductor.Nturbo; // Meter} {Länge des Spulendrahtes, Turbo-Spule
            RI = param.General.rho * (DLI) / (Math.PI / 4.0 * param.Inductor.DD * param.Inductor.DD); // Ohm} {Ohm`scher Widerstand des Spulendrahtes, Input-Spule
            RT = param.General.rho * (DLT) / (Math.PI / 4.0 * param.Inductor.DD * param.Inductor.DD); // Ohm} {Ohm`scher Widerstand des Spulendrahtes, Turbo-Spule
            BreiteI = param.Inductor.Nebeninput * param.Inductor.DD; HoeheI = param.Inductor.Ueberinput * param.Inductor.DD; // Breite und Höhe des Input-Spulenlörpers
            BreiteT = param.Inductor.Nebenturbo * param.Inductor.DD; HoeheT = param.Inductor.Ueberturbo * param.Inductor.DD; // Breite und Höhe des Turbo-Spulenlörpers
            fkI = Math.Sqrt(HoeheI * HoeheI + 4.0 / Math.PI * 2.0 * yo * 2.0 * zo) / HoeheI; // Korrekturfaktor zur Induktivität der kurzen Inpuz-Spule nach *1 von S.16
            fkT = Math.Sqrt(HoeheT * HoeheT + 4.0 / Math.PI * 2.0 * yo * 2.0 * zo) / HoeheT; // Korrekturfaktor zur Induktivität der kurzen Turbo-Spule nach *1 von S.16
            PrintLine("Induktivitäts-Korrektur: fkI = " + fkI + ",  fkT = " + fkT);
            LI = Physics.Muo * (2.0 * yo + BreiteI) * (2.0 * zo + BreiteI) * param.Inductor.Ninput * param.Inductor.Ninput / (HoeheI * fkI); // Geometrische Mittelung => Induktivität der Input-Spule
            LT = Physics.Muo * (2.0 * yo + BreiteT) * (2.0 * zo + BreiteT) * param.Inductor.Nturbo * param.Inductor.Nturbo / (HoeheT * fkT); // Geometrische Mittelung => Induktivität der Turbo-Spule
            omT = 1.0 / Math.Sqrt(LT * param.General.CT);  // Kreis-Eigenfrequenz des Turbo-Spulen-Schwingkreises aus LT & CT
            TT = 2.0 * Math.PI / omT;        // Schwingungsdauer des Turbo-Spulen-Schwingkreises aus LT & CT.
            Mmag = param.General.rhoMag * (Math.PI * param.PermanentMagnet.MEro * param.PermanentMagnet.MEro) * (2.0 * param.PermanentMagnet.MEyo);// Masse des Magneten}                     {Rotation des Magneten um Querachse !!
            J = Mmag / 4.0 * (param.PermanentMagnet.MEro * param.PermanentMagnet.MEro + 4.0 * param.PermanentMagnet.MEyo * param.PermanentMagnet.MEyo / 3.0); // Trägheitsmoment des Magneten bei Rotation, siehe *2 von S.13 und Dubbel S.B-32
            //  Sonstige, auch abgeleitete (aus den obigen Parametern berechnete) Größen:
            omAn = param.General.UmAn / 60.0 * 2.0 * Math.PI;  // Rotierender Magnet: Winkelgeschwindigkeit (rad/sec.), Startdrehzahl
            UmSec = param.General.UmAn / 60.0;      // Rotierender Magnet: Umdrehungen pro Sekunde, Startdrehzahl

            //  Anzeige der Werte:
            PrintLine("*******************************************************************************");
            PrintLine("Anzeige einiger auszurechnender Parameter:");
            PrintLine("Länge des Spulendrahtes, Input-Spule: " + DLI + " m");
            PrintLine("Länge des Spulendrahtes, Turbo-Spule: " + DLT + " m");
            PrintLine("Ohm`scher Widerstand des Input-Spulendrahts: RI = " + RI + " Ohm");
            PrintLine("Ohm`scher Widerstand des Turbo-Spulendrahts: RT = " + RT + " Ohm");
            PrintLine("Induktivität der Input-Spule, ca.: LI = " + LI + " Henry");
            PrintLine("Induktivität der Turbo-Spule, ca.: LT = " + LT + " Henry");
            PrintLine("Eigen-Kreisfreq des Turbo LT-CT-Schwinkreises: omT = " + omT + " Hz (omega)");
            PrintLine("=> Schwingungsdauer TT = 2*Math.PI/omT = " + TT + "sec.");
            PrintLine("Magnet: Start-Winkelgeschw.: omAn = " + omAn + " rad/sec");
            PrintLine("Magnet: Startdrehzahl, Umdr./sec.: UmSec = " + UmSec + " Hz");
            PrintLine("Masse des Magnet = " + Mmag + " kg");
            PrintLine("Trägheitsmoment Magnet bei QUER-Rotation" + J + " kg*m^2");
            PrintLine("Gesamtdauer der Betrachtung: " + param.AnzP * param.dt + " sec.");
            PrintLine("Excel-Export: " + param.PlotAnfang * param.dt + "..." + param.PlotEnde * param.dt + " sec., Step " + param.PlotStep * param.dt + " sec.");
            PrintLine("Das sind " + (param.PlotEnde - param.PlotAnfang) / param.PlotStep + " Datenpunkte (also Zeilen).");
            if (((param.PlotEnde - param.PlotAnfang) / param.PlotStep) > AnzPmax)
            {
                PrintLine(); PrintLine("FEHLER: Zu viele Datenpunkte.");
                PrintLine("So viele Datenpunkte können in Excel nicht dargestellt werden.");
                PrintLine("=> Berechnung wird hier GESTOPPT.");
                Environment.Exit(1);
            }

            //  Wait; 
            //  Hilfsarbeiten: Für die Vorbereitungen brauche ihch AnzP=360, danach brauche ich wieder den eingegebenen Wert.
            ulong AnzPmerk = param.AnzP; // Merken des Wertes für später
            param.AnzP = 360;      // Eine Umdrehung in Winkel-Grad-Schritten

            //  Hier beginnt das Rechenprogramm.
            // Teil 1: Eine Test-Berechnung der Drehmoment-Wirkung der beiden Spulen auf den Magneten:
            PrintLine();  // Wir beginnen mit der Bestimmung des Magnetfeld des beiden Spulen an einem beliebigen Aufpunkt
            PrintLine("Hier steht vorerst zu Testzwecken die Feldberechnung der Input- und Turbo-Spule");
            this.MagfeldInputBerechnen(param, -0.00, 0.01, 0.01, 1.0);  // drei kartesische Komponenten für Aufpunkt und Strom = 1.0 Ampere
            PrintLine("B_Input_x,y,z:" + BIx + ", " + BIy + ", " + BIz + " T");
            this.MagfeldTurboBerechnen(param, +0.00, 0.01, 0.01, 1.0);  // drei kartesische Komponenten für Aufpunkt und Strom = 1.0 Ampere
            PrintLine("B_Turbo_x,y,z:" + BTx + ", " + BTy + ", " + BTz + " T");
            merk = Math.Sqrt((2 * yo * Spsw * 2 * zo * Spsw) + Physics.Sqr(xo * Spsw)); merk = merk * merk * merk;
            PrintLine("Vgl->Input: Runde Leiterschleife, Feld im Ursprung: " + Physics.Muo * param.Inductor.Ninput * 1.0 * 2 * yo * Spsw * 2 * zo * Spsw / 2 / merk + " T");
            PrintLine("Vgl->Turbo: Runde Leiterschleife, Feld im Ursprung: " + Physics.Muo * param.Inductor.Nturbo * 1.0 * 2 * yo * Spsw * 2 * zo * Spsw / 2 / merk + " T");
            // Die Berechnung des Magnetfeldes der beiden Spulen (Input & Turbo) ist jetzt getestet und funktioniert.
            if (!schonda)
            {         // Das ist nur eine Kontrolle.
                for (int i = 0; i < 360; i++)
                {
                    simPrep[i].VoltageInducedInput = Drehmoment(param, 1.0 * i / 180.0 * Math.PI);
                    PrintLine(i + "Grad => Drehmoment-Komponente: Mz = " + simPrep[i].VoltageInducedInput); // Das Argument ist der Winkel der Magnetstellung "alpha"
                }

                ExcelAusgabe(param, simPrep, "drehmom.dat", 14);   // Optionales Upgm. zur Datenausgabe nach Excel.
                PrintLine("Damit ist die Drehmoments-Berechnung des Magneten geschafft.");
            }

            // Teil 2: Ausprobieren der Flußberechnung durch die beiden Spulen unter Magnet-Drehung (führt später zur induz. Spannung):
            if (!schonda)
            {
                PrintLine("Es folgt die Berechnung des magnetischen Flusses für Geometrie \"03\"");
                this.MagnetDrehen(00); // Drehwinkel in Grad angeben, 0...360
                this.GesamtflussDurchInputSpule();
                PrintLine("Gesamtfluss durch Input-Spule: " + PsiGES + " T*m^2");
                this.MagnetDrehen(01); // Drehwinkel in Grad angeben, 0...360
                this.GesamtflussDurchInputSpule();

                PrintLine("Gesamtfluss durch Input-Spule: " + PsiGES + " T*m^2");
                PrintLine("-----------------------");
                this.MagnetDrehen(00); // Drehwinkel in Grad angeben, 0...360
                this.GesamtflussDurchTurboSpule();
                PrintLine("Gesamtfluss durch Turbo-Spule: " + PsiGES + " T*m^2");
                this.MagnetDrehen(01); // Drehwinkel in Grad angeben, 0...360
                this.GesamtflussDurchTurboSpule();
                PrintLine("Gesamtfluss durch Turbo-Spule: " + PsiGES + " T*m^2");
                PrintLine("-----------------------");
            }
            // Ergebnis bis hier: Die Differenz zwischen beiden erlaubt die Berechnung der induzierten Spannung

            //  Test: Einmal den Magneten drehen und den magnetischen Fluss / die induzierte Spannung messen:
            // Zum Test stehen 360 Zeiteinheiten = 360*dt = 36 Millisekunden für eine Umdrehung, entsprechend 1666.666 U/min
            if (!schonda)
            {
                PrintLine("Zuerst die Input-Spule:");

                for (int i = 0; i < 360; i++)  // Zuerst probiere ich's mit der Input-Spule
                {
                    if (this.Cancel)
                    {
                        return CreateEmptySimulationResults();
                    }

                    simPrep[i].phi = i; // Angabe in Grad
                    this.MagnetDrehen(simPrep[i].phi);

                    this.GesamtflussDurchInputSpule();  // setzt auf "PsiGES" den Ergebnis-Wert ab.

                    simPrep[i].PSIinput = PsiGES; // Dies ist der magnetische Fluß durch die Input-Spule
                    PrintLine("phi = " + simPrep[i].phi + " grad  =>  magn. ges. Fluss = " + simPrep[i].PSIinput + " T*m^2");

                    if (i == 0)
                    {
                        simPrep[i].UindInput = 0;
                    }

                    if (i > 0)
                    {
                        simPrep[i].UindInput = -param.Inductor.Ninput * (simPrep[i].PSIinput - simPrep[i - 1].PSIinput) / param.dt;
                    }

                    simPrep[i].KG = 0;
                    simPrep[i].KH = simPrep[i].PSIinput;
                    simPrep[i].KI = simPrep[i].UindInput;  // Zur Excel-Ausgabe weiterleiten

                }

                PrintLine("---------------------------------");
                PrintLine("Danach die Turbo-Spule:");

                for (int i = 0; i < 360; i++)  // Danach probiere ich's auch noch mit der Turbo-Spule
                {
                    if (this.Cancel)
                    {
                        return CreateEmptySimulationResults();
                    }

                    simPrep[i].phi = i; // Angabe in Grad

                    this.MagnetDrehen(simPrep[i].phi);
                    this.GesamtflussDurchTurboSpule();  // setzt auf "PsiGES" den Ergebnis-Wert ab.

                    simPrep[i].PSIturbo = PsiGES; // Dies ist der magnetische Fluß durch die Turbo-Spule

                    PrintLine(String.Format("phi = {0} grad  =>  magn. ges. Fluss = {1} T*m^2", simPrep[i].phi, simPrep[i].PSIturbo));

                    if (i == 0)
                    {
                        simPrep[i].UindTurbo = 0;
                    }

                    if (i > 0)
                    {
                        simPrep[i].UindTurbo = -param.Inductor.Nturbo * (simPrep[i].PSIturbo - simPrep[i - 1].PSIturbo) / param.dt;
                    }

                    simPrep[i].KJ = 0;
                    simPrep[i].MagneticFluxTurbo = simPrep[i].PSIturbo;
                    simPrep[i].MageneticFluxInput = simPrep[i].UindTurbo;  // Zur Excel-Ausgabe weiterleiten
                    simPrep[i].VoltageInducedTurbo = 0;
                    simPrep[i].VoltageInducedInput = simPrep[i].VoltageInducedInput; // Am Ende noch zwei Leerspalten
                }

                // Jetzt muß man das Signal noch gegen Rauschen glätten:
                // FourierDatenspeicherung(sim, x => x.PSIturbo);
                AdvancedMath.FourierEntwicklung(param, simPrep, x => x.PSIturbo, out B1T, out B2T, out B3T, out B4T, out B5T);

                // FourierDatenspeicherung(sim, x => x.PSIinput);
                AdvancedMath.FourierEntwicklung(param, simPrep, x => x.PSIinput, out B1I, out B2I, out B3I, out B4I, out B4I);

                // Kontroll-Output der geglätteten Fluß-Werte in den ersten beiden Kolumnen des Excel-Datensatzes:
                for (int i = 0; i < 360; i++)
                {       // FlussI und FlussT gibt den geglätteten magnetischen Fluß durch die Spulen an.
                    simPrep[i].KJ = FlussI(1.0 * i / 360.0 * 2.0 * Math.PI);     // Der Lagewinkel des Magneten wird in Radianten angegeben.
                    simPrep[i].VoltageInducedTurbo = FlussT(1.0 * i / 360.0 * 2.0 * Math.PI);     // Der Lagewinkel des Magneten wird in Radianten angegeben.
                }
            }

            // Die Drehmoments-Berechnung absorbiert noch zu viel CPU-Zeit, um mit sehr feiner Zeit-Schrittweite rechnen zu können.
            // Daher entwickele ich jetzt auch eine Fourier-Reihe zur Beschleunigung der Drehmoments-Berechnung:
            if (!schonda)
            {
                qpoT = 1; qpoI = 0; // Schnell-Berechnungs-Kalibrierung nur für Turbo-Spule, 1A, aber ohne Input-Spule
                PrintLine("Drehmoment in einen Sinus-Term umrechnen, zur späteren CPU-Zeit Ersparnis:");
                for (int i = 0; i < 360; i++)
                {       // Das gesamt-Drehmoment, das der Magnet im Feld der beiden Spulen (Input&Turbo) aufnimmt.
                    simPrep[i].VoltageInducedInput = Drehmoment(param, 1.0 * i * 2.0 * Math.PI / 360.0); // Der Lagewinkel des Magneten wird in Radianten angegeben.
                    //Print(".");  // PrintLine(KN[i]);
                }

                AdvancedMath.SinusEntwicklung_fuer_Drehmoment(param, simPrep, x => x.VoltageInducedInput, out B1dreh, out phase);

                PrintLine("Drehmom-Ampl: " + B1dreh + " und Phase: " + phase);
                // Kontrolle, ob die Schnell-Drehmoments-Berechnung richtige Ergebnisse liefert:
                for (int i = 0; i < 360; i++)
                {
                    simPrep[i].KG = this.SchnellDrehmoment(1.0 * i * 2.0 * Math.PI / 360.0); // Der Lagewinkel des Magneten wird in Radianten angegeben.
                }
            }

            // Daten abspeichern, falls eine Parameter-Konfiguration vorliegt:
            if (!schonda && !this.Cancel)
            {
                this.DatenSpeichern(param);
            }

            // Damit ist die Vorbereitung beendet.

            // Ich kontrolliere jetzt, ob alle Parameter und Daten mit und ohne "schonda" angekommen sind:
            for (int i = 0; i < 360; i++)
            {       // FlussI und FlussT gibt den geglätteten magnetischen Fluß durch die Spulen an.
                simPrep[i].KJ = FlussI(1.0 * i * 2.0 * Math.PI / 360.0); // Magnetischer Fluss durch Input-Spule, Winkel des Magneten in Radianten
                simPrep[i].VoltageInducedTurbo = FlussT(1.0 * i * 2.0 * Math.PI / 360.0); // Magnetischer Fluss durch Turbo-Spule, Winkel des Magneten in Radianten
            }

            for (int i = 0; i < 360; i++)
            {
                simPrep[i].KG = this.SchnellDrehmoment(1.0 * i * 2.0 * Math.PI / 360.0); // Drehmoment auf den Magneten, Winkel des Magneten in Radianten
            }

            ExcelAusgabe(param, simPrep, "kontroll.dat", 14);  // Optionales Upgm. zur Datenausgabe nach Excel.
            // Hilfsarbeiten: Zum Lösen der Dgl. brauche ich wieder vorgegebene Anzahl von Iterationsschritten:
            param.AnzP = AnzPmerk;
            // Damit stehen jetzt alle Daten zur DFEM-Berechnung bereit.
            PrintLine("*******************************************************************************");


            var sim = new SimulationDataItem[AnzPmax];

            for (int i = 0; i < sim.Length; i++)
            {
                sim[i] = new SimulationDataItem();
            }

            // Initialisierung für die Suche der Maximalwerte zur Strom-, Drehzahl- und Spannungsangabe bei der Auslegung:
            QTmax = 0; QImax = 0; QpTmax = 0; QpImax = 0; QppTmax = 0; QppImax = 0; phipomax = 0;
            Wentnommen = 0; // Initialisierung der entnommenen Energie am Lastwiderstand
            Ereib = 0;      // Initialisierung der entnommenen mechanischen Energie über Reibung

            // Initialisierung der Referenz für das Input-Spannungs-Signal:
            steigtM = false; steigtO = false;
            Ezuf = 0;  // Initialisierung der zugeführten Energie über die Input-Spannung
            LPP = 0;   // Initialisierung der Anzahl der Datenpunkt für den Excel-Plot

            //  Damit ist alles vorbereitet, und ich kann jetzt anfangen, in die Lösung der System-Differentialgleichung zu gehen:
            //  Also kommt jetzt der Rechen-Kern:
            //  Zuerst: Anfangsbedingungen einsetzen:
            phio = 0; phipo = omAn;  // phippo:=0; //Anfangsbedingungen der Magnetrotation (mechanisch)
            // Starten mit gegebener Anfangs-Winkelgeschwindigkeit ist vorgesehen.
            qoT = param.General.CT * param.General.Uc;
            qpoT = param.General.Il;     // qppoT:=0; //Anfangsbedingungen Turbo-Spule (elektrisch)
            // Ein geladener Kondensator im Schwingkreis der Turbo-Spule ist vorgesehen.
            qoI = 0; qpoI = 0; qppoI = 0;   // Anfangsbedingungen Input-Spule (elektrisch), zunächst inaktiv
            //  Beim "nullten" Schritt ist der alten Schritt nicht der "minus erste", sondern auch der nullte Schritt:
            // phim:=phio;}{phipm:=phipm;} {phippm:=phippm;
            qmT = qoT;    // qpmT:=qpmT;}   {qppmT:=qppmT;
            // qmI:=qoI;}  {qpmI:=qpmI;}   {qppmI:=qppmI;

            this.OnInitializeProgressReport((param.PlotEnde - param.PlotAnfang) / param.PlotStep);

            //  Jetzt das eigentliche Lösen der Differential-Gleichung:}    {Zunächst noch ohne Input-Spule !!
            for (ulong i = 0; i < param.AnzP; i++)
            {
                if (this.Cancel)
                {
                    break;
                }

                // Initialisierung der Referenz für das Input-Spannungs-Signal:
                if (i == 0)
                {
                    Iumk = 0;
                }

                if (i >= 1) // Input-Spannungs-Referenz hier an der Turbo-Spule festmachen.
                {
                    steigtM = steigtO;  // alten Flanken-Steigungs-Zustand merken
                    if (qoT > qmT) steigtO = true;
                    if (qoT < qmT) steigtO = false;
                    if (steigtM && !steigtO) Iumk = i;
                }

                // Aktueller Moment der Analyse, laufende Zeit in Sekunden, "Jetzt-Schritt":
                Tjetzt = i * param.dt;
                // Für den neuen Schritt, wird der alte Vorgänger-Schritt zum vorletzten Schritt heruntergezählt:
                phim = phio; phipm = phipo; /* phippm:=phippo;*/   // Drehbewegung
                qmT = qoT; qpmT = qpoT;     /* qppmT:=qppoT;*/     // Turbo-Spule
                qmI = qoI; qpmI = qpoI; qppmI = qppoI;       // Input-Spule
                // Und jetzt rechne ich den neuen Schritt aus:
                // Zuerst die Drehung des Magneten, Drehmoment kommt aus Spulenströmen. Vorhandene Drehmoments-Berechnung benutzen:
                /* KK */
                phippo = this.SchnellDrehmoment(phim) * qpoT / J;  // Upgm "Schnell_Drehmoment" ist skaliert ITurbo=1A & IInput=0A, geht linear mit dem Turbo-Strom.
                // !! Alle Zeilen mit )ppelten Ausrufezeichen dienen der zusätzlichen Aufnahme der Input-Spule.
                // !! phippo:=Drehmoment(phim)/J;}               {Vollwertige Drehmoment-Berechnung mit Turbo-Spule und Input-Spule, geht langsam.
                // Für phippo muß ich eine der beiden vorangehenden Alternativen verwenden, je nach dem ob nur die Turbo-Spule aktiv ist, oder auch die Input-Spule.
                // Falls die Input-Spule auch aktiv ist, dann soll ich "schonda" von Anfang an auf "false" setzen und immer die komplette Vorbereitung durchrechnen.
                // Alle mit "GG" kommentierten Zeilen dienen einer geschwindigkeits-proportionalen mechanischen Leistungsentnahme:
                /* GG */
                if (i == 1) cr = param.General.crAnfang;            // Geschwindigkeits-proportionaler Reibungs-Koeffizient
                /* GG */
                if (i > 1) cr = this.ReibungNachregeln(param.General);  // Dieser wird nachgeregelt, um eine konstante Drehzahl einzustellen für stabilen Betrieb der Maschine.
                /* GG */
                if (phipo > 0) phippo = phippo - cr * phipm / J; //Die neg. Beschleunigung wirkt immer der Geschwindigkeit entgegen.
                /* GGif( phipo=0 ) phippo:=phippo;
                /* GG */
                if (phipo < 0) phippo = phippo + cr * phipm / J; // Die neg. Beschleunigung wirkt immer der Geschwindigkeit entgegen.
                /* GG */
                // Jetzt ist die geschwindigkeits-proportionale Reibung berechnet.
                //if ((i % 100000) == 0) Print(".");
                phipo = phipm + phippo * param.dt;                   // 1. Integrationsschritt, ohne mechanische Reibung
                phio = phim + phipo * param.dt;                      // 2. Integrationsschritt
                /* GG */
                Preib = cr * phipm * phipo;       // Leistungsentnahme zur geschwindigkeits-proportionalen Reibung
                /* GG */
                Ereib = Ereib + Preib * param.dt;       // Gesamte über Reibung entnommene Leistung
                // Dann die Turbo-Spule. Gedämpfte elektrische Schwingung, dazu induzierte Spannung aufgrund Magnet-Drehung:
                /* FF */
                qppoT = -1.0 / (LT * param.General.CT) * qmT - (RT + param.General.Rlast) / LT * qpmT; //Dgl. der gedämpfte Schwingung.
                UinduzT = -param.Inductor.Nturbo * (FlussT(phio) - FlussT(phim)) / param.dt;// Wirkung durch die induzierte Spannung (aufgrund der Magnet-Drehung) hinzunehmen
                qppoT = qppoT - UinduzT / LT;                  // Wirkung der induzierte Spannung auf die zweite Ableitung von q, also "qppoT")
                /* ??*/
                qpoT = qpmT + qppoT * param.dt; /*-Rlast/(2*LT)*qpmT*dt;*/ // 1. Integrationsschritt, nach *5 von S.6 im alten Skript bzw *1 von S.14 im neunen Skript
                qoT = qmT + qpoT * param.dt;                         // 2. Integrationsschritt, nach *3 & *4 von S.6 im alten Skript bzw *1 von S.14 im neunen Skript
                /* Dann die Input-Spule:  */
                UinduzI = 0;
                qoI = qmI; qpoI = qpmI; qppoI = qppmI;          // Die Input-Spule macht noch gar nichts, sie sMath.PIelt jetzt noch nicht mit.
                // Und wenn die Input-Spule noch mitspielt, muß ich die nachfolgenden fünf Zeilen zur Input-Spule aktivieren:
                // !! qppoI:=-1/(LI*CI)*qmI-RI/LI*qpmI+U7/LI;}      {Dgl. der gedämpfte Schwingung, dazu Störfunktion für Input-Spannung in den Input-Schwingkreis
                // !! UinduzI:=-Ninput*(FlussI(phio)-FlussI(phim))/dt;} {Wirkung durch die induzierte Spannung (aufgrund der Magnet-Drehung) hinzunehmen
                // !! qppoI:=qppoI-UinduzI/LI;}                     {Wirkung der induzierte Spannung auf die zweite Ableitung von q, also "qppoT")
                // !! qpoI:=qpmI+qppoI*dt;}                         {1. Integrationsschritt, nach *5 von S.6 im alten Skript bzw *1 von S.14 im neunen Skript
                // !! qoI:=qmI+qpoI*dt;}                            {2. Integrationsschritt, nach *3 & *4 von S.6 im alten Skript bzw *1 von S.14 im neunen Skript
                Pzuf = U7(i, Iumk);// *qpoI}                             {Über die Input-Spannung zugeführte Leistung
                Ezuf = Ezuf + Pzuf * param.dt;                          // Zugeführte Energie  über die Input-Spannung
                // Achtung: Die Drehmoments-Schnell-Berechnung "phippo" geht so noch nicht für Turbo-Input-Spule. Dazu muß ich noch die Ströme in die einzelnen Upgme durchreichen.

                // Jetzt muß ich noch die Maximalwerte für Strom, Spannung und Drehzahl der Auslegung bestimmen:
                if (Math.Abs(qoT) > QTmax) QTmax = Math.Abs(qoT);           // Maximum der Ladung im Turbo-Kondensator festhalten
                if (Math.Abs(qoI) > QImax) QImax = Math.Abs(qoI);           // Maximum der Ladung im Input-Kondensator festhalten
                if (Math.Abs(qpoT) > QpTmax) QpTmax = Math.Abs(qpoT);       // Maximum des Stroms in der Turbo-Spule festhalten
                if (Math.Abs(qpoI) > QpImax) QpImax = Math.Abs(qpoI);       // Maximum des Stroms in der Input-Spule festhalten
                if (Math.Abs(qppoT) > QppTmax) QppTmax = Math.Abs(qppoT);   // Maximum des Ipunkt in der Turbo-Spule festhalten
                if (Math.Abs(qppoI) > QppImax) QppImax = Math.Abs(qppoI);   // Maximum des Ipunkt in der Input-Spule festhalten
                if (Math.Abs(phipo) > phipomax) phipomax = Math.Abs(phipo); // Maximum der Winkelgeschwindigkeit des Magneten festhalten
                Wentnommen += param.General.Rlast * qpoT * qpoT * param.dt;  // Summierung der am Lastwiderstand im Turbo-Schwingkreis entnommenen Gesamtenergie

                // Ggf. muß ich jetzt einen Plot-Punkt für's Excel ausgeben:
                if ((i >= param.PlotAnfang) && (i <= param.PlotEnde)) // diese Punkte ins Excel plotten
                {
                    if (((i - param.PlotAnfang) % (param.PlotStep)) == 0)
                    {
                        znr = (i - param.PlotAnfang) / param.PlotStep;

                        sim[znr].Zeit = Tjetzt;                                                         // Zeitpunkt für Excel abspeichern.
                        sim[znr].Q = qoT;                                                               // Turbo-Spule, im Array (und nur )rt) ohne Index "T".
                        sim[znr].Qp = qpoT;
                        sim[znr].Qpp = qppoT;
                        sim[znr].QI = qoI;                                                              // Input-Spule
                        sim[znr].QpI = qpoI;
                        sim[znr].QppI = qppoI;
                        sim[znr].phi = phio / AdvancedMath.TwoPI; /* Umdrehungen */                    // Drehung des Magneten
                        sim[znr].phip = phipo / AdvancedMath.TwoPI * 60; /* U/min */
                        sim[znr].phipp = phippo / AdvancedMath.TwoPI * 3600; /*U/min/min*/
                        sim[znr].MagneticFluxTurbo = FlussT(phio);                                     // Magnetischer Fluß durch die Spulen
                        sim[znr].MageneticFluxInput = FlussI(phio);
                        sim[znr].VoltageInducedTurbo = UinduzT;
                        sim[znr].VoltageInducedInput = UinduzI;                                         // In den Spulen induzierte Spannung
                        sim[znr].EnergyInputCoil = 1.0 / 2.0 * LT * qpoT * qpoT;                        // Energie in der Input-Spule
                        sim[znr].EnergyTurboCoil = 1.0 / 2.0 * LI * qpoI * qpoI;                        // Energie in der Turbo-Spule
                        sim[znr].EnergyInputCapacitor = 1.0 / 2.0 * qoT * qoT / param.General.CT;       // Energie im Input-Kondensator
                        sim[znr].EnergyTurboCapacitor = 1.0 / 2.0 * qoI * qoI / param.General.CI;       // Energie im Turbo-Kondensator
                        sim[znr].EnergyMagneticRotation = 1.0 / 2.0 * J * phipo * phipo;                // Energie der Magnet-Rotation
                        sim[znr].EnergySystem = sim[znr].EnergyInputCoil + sim[znr].EnergyTurboCoil + sim[znr].EnergyInputCapacitor + sim[znr].EnergyTurboCapacitor + sim[znr].EnergyMagneticRotation;     // Gesamt-Energie im System
                        sim[znr].PowerElectric = param.General.Rlast * qpoT * qpoT;                                // Am Lastwiderstand entnommene Leistung, nur Turbo-seitig
                        sim[znr].VoltageInputCoil = U7(i, Iumk);                                                      // Kontrolle der Input-Spannung im Input-Schwingkreis
                        sim[znr].PowerAddedByInput = Pzuf;                                        // Zugeführte Leistung über die Input-Spannung
                        sim[znr].FrictionCoefficient = cr;                                          // geregelter Reibungskoeffizient zur mechanischen Leistungsentnahme
                        sim[znr].MechanicalFrictionPowerOut = Preib;                                       // Entnommene mechanische Leistung, emuliert durch geschwindigkeits-proportionale Reibung
                        LPP = znr;                                             // Letzter Plot-Punkt; Wert wird für Datenausgabe benutzt -> ExcelLangAusgabe

                        if (znr % this.ProgressReportStepCount == 0)
                        {
                            this.ProgressReport(znr, sim[znr]);
                        }
                    }
                }

                AnfEnergie = sim[0].PowerElectric;         // Anfangs-Gesamt-Energie im System
                EndEnergie = sim[LPP].PowerElectric;       // End-Gesamt-Energie im System
            }

            if (this.Cancel)
            {
                PrintLine("");
                PrintLine("Simulation wurde abgebrochen");

                return CreateEmptySimulationResults(sim);
            }

            var summary = this.CreateSimulationResultSummary(param, sim);

            this.ExcelLangAusgabe(param, sim, "test.dat", 25);

            var result = new SimulationResult();
            result.Summary = summary;
            result.SimulationData = sim;

            PrintResult(result.Summary, PrintLine);

            return result;
        }

        private static SimulationResult CreateEmptySimulationResults(SimulationDataItem[] simData = null)
        {
            if (simData == null)
            {
                simData = new SimulationDataItem[0];
            }

            return new SimulationResult
                {
                    SimulationData = simData,
                    Summary = new SimulationResultSummary()
                };
        }

        private SimulationResultSummary CreateSimulationResultSummary(SimulationParameters param, SimulationDataItem[] sim)
        {
            var summary = new SimulationResultSummary();

            summary.LPP = LPP;
            summary.UmaxI = QImax / param.General.CI;
            summary.UmaxT = QTmax / param.General.CT;
            summary.ImaxI = QpImax;
            summary.ImaxT = QpTmax;
            summary.UmaxSI = LI * QppImax;
            summary.UmaxST = LT * QppTmax;
            summary.DrehzahlMaxRad = phipomax;
            summary.DrehzahlMaxUmin = phipomax / 2 / Math.PI * 60;
            summary.DrehzahlEndeUmin = sim[LPP].phip / 2 / Math.PI * 60;
            summary.AnfEnergie = AnfEnergie;
            summary.EndEnergie = EndEnergie;
            summary.EnergieZunahme = (EndEnergie - AnfEnergie);
            summary.LeistungsAenderung = (EndEnergie - AnfEnergie) / (param.AnzP * param.dt);
            summary.PentnommenRlast = Wentnommen;
            summary.PentnommenRlastMittel = Wentnommen / (param.AnzP * param.dt);
            summary.Ezugefuehrt = Ezuf;
            summary.Pzugefuehrt = Ezuf / (param.AnzP * param.dt);
            summary.Emech = Ereib;
            summary.Pmech = Ereib / (param.AnzP * param.dt);
            summary.DauerSimulation = (param.AnzP * param.dt);

            return summary;
        }

        private static void PrintResult(SimulationResultSummary summary, Action<string> printLine)
        {
            printLine("");
            printLine("Anzahl Datenpunkte für Excel-Plot: LPP = " + summary.LPP);
            printLine("");
            printLine("Anzeigen der Amplituden der Auslegung: (nicht Effektivwerte, sondern Spitze)");
            printLine("Input-Kondensator, Spannung, UmaxI = " + summary.UmaxI + " Volt");   // Maximum der Ladung im Input-Kondensator
            printLine("Turbo-Kondensator, Spannung, UmaxT = " + summary.UmaxT + " Volt");   // Maximum der Ladung im Turbo-Kondensator
            printLine("Input-Schwingkreis,   Strom, ImaxI = " + summary.ImaxI + " Ampere");   // Maximum des Stroms in der Input-Spule
            printLine("Turbo-Schwingkreis,   Strom, ImaxT = " + summary.ImaxT + " Ampere");   // Maximum des Stroms in der Turbo-Spule
            printLine("Input-Spule,       Spannung, UmaxSI= " + summary.UmaxSI + " Volt"); // Maximum des Ipunkt in der Input-Spule
            printLine("Turbo-Spule,       Spannung, UmaxST= " + summary.UmaxST + " Volt"); // Maximum des Ipunkt in der Turbo-Spule
            printLine("Maximale Magnet-Rotationsdrehzahl  = " + summary.DrehzahlMaxRad + " rad/sec");// Maximum der Winkelgeschwindigkeit des Magneten
            printLine("Maximale Magnet-Rotationsdrehzahl= " + summary.DrehzahlMaxUmin + " U/min."); // Maximum der Winkelgeschwindigkeit des Magneten
            printLine("Am Ende erreichte End-Drehzahl   = " + summary.DrehzahlEndeUmin + " U/min.");
            printLine("");
            printLine("Anfangs-Energie im System:    " + summary.AnfEnergie + " Joule");
            printLine("End-Energie im System:        " + summary.EndEnergie + " Joule");
            printLine("Energie-Zunahme im System:    " + summary.EnergieZunahme + " Joule");
            printLine("Leistungs-Änderung im System: " + summary.LeistungsAenderung + " Watt");
            printLine("Am Lastwiderstand entnommene Gesamtenergie = " + summary.PentnommenRlast + " Joule");
            printLine("entsprechend einer mittleren entnommenen Leistg: " + summary.PentnommenRlastMittel + " Watt");
            printLine("Über Input-Spannung zugeführte Gesamt-Energie: " + summary.Ezugefuehrt + " Joule");
            printLine("entsprechend einer mittleren zugeführten Leistg: " + summary.Pzugefuehrt + " Watt");
            printLine("Gesamte mechanisch entnommene Energie = " + summary.Emech + " Joule");
            printLine("entsprechend einer mittleren Leistung = " + summary.Pmech + " Watt");
            printLine("bei einer Betrachtungs-Dauer von " + summary.DauerSimulation + " sec.");
            printLine("");
            printLine("Fertig gerechnet -> Adele.");
        }

        private void Print(string message)
        {
            if (this.LogCallback != null)
            {
                this.LogCallback(message);
            }
        }

        private void PrintLine()
        {
            if (this.LogCallback != null)
            {
                this.LogCallback(Environment.NewLine);
            }
        }

        private void PrintLine(string message)
        {
            if (this.LogCallback != null)
            {
                this.LogCallback(message + Environment.NewLine);
            }
        }

        private void OnInitializeProgressReport(ulong stepCount)
        {
            if (this.InitializeProgressReport != null)
            {
                this.InitializeProgressReport(stepCount);
            }
        }

        public Action<string> LogCallback { get; set; }
        public Action<ulong> InitializeProgressReport { get; set; }
        public Action<ulong, SimulationDataItem> ProgressReport { get; set; }
        public ulong ProgressReportStepCount { get; set; }
    }
}
